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1. Introduction 

This is the final report for the period of performance 7/1/88-6/30/92. The project 
was originally proposed for three years. Funds were not made available for the third year. 
A no-cost extension of the second year was requested in order to complete as much of the 
original project as possible. 

2. Summary 

The following tasks were described in the original and continuation proposals. 


Task 1: Order of Magnitude Estimates 

Task 2 : Thermo-capillary Convection - Two Dimensional (fixed planar surface) 

Task 3: Thermo-capillary Convection - Three Dimensional and Axisymmetric 

Task 4: Liquid Bridge/ Floating-Zone Sensitivity 

Task 6: Transport in Closed Containers 

Task 7 : Interaction: Design and Development Stages 

Task 8: Interaction: Testing Flight Hardware 

Task 9: Reporting 


Work was completed on all aspects of these tasks despite the lack of funding for the 
third year. The only part of the origin^ proposal which was not addressed was the effect of 
vibration on diffusion experiments. 

A review of work related to experiment sensitivity to residual acceleration has been 
published and a copy of the article is attached in Appendix 1. In the third semi-annual 
report we stated that Task 2 should complete by June-July of 1990 (i.e. earlier than 
originally anticipated). Since that report was issued NASA MSFC extended the P.O.P. of 
the second year’s work to June 1991 as continuation funds are not yet available. We 
completed Task 2 during this period and a paper describing the results is in preparation for 
journal publication. (Task 2 involves an analysis of the residual acceleration sensitivity of 
thermocapillary driven flow in a 2-D rectangular region. The free surface is held planar. ) 

The work is summarized in Section 4. 

Task 3, which involves the examination of axisymmetric thermocapillary 
convection has been started early in the process of developing a method for the examination 
of non-isothermal liquid bridge. A code for the axisymmetric isothermal and non- 
isothermal (i.e. thermocapillary driven convection is included along with a deformable free 
surface) liquid bridge sensitivity analysis has been developed and a comparison with 1-D 
results has been made. These results will be reported in Section 3. 

Five papers, and one extended abstract, summarizing work carried out entirely (or 
partly*) under the VIT project have been accepted for publication: 

1. Microgravity Experiment Sensitivity to Residual Accelerations: A Review, 
Microgravity Science and Technology ,m, No.2, 52-68, 1990. (See Appendix 1) 

2. The Sensitivity of a Liquid Bridge to Axial Vibration, with Y. Q. Zhang, 
Physics of Fluids, A 2, 1966-1974, 1990. (See Appendix 2) 


3.* Bridgman Crystal Growth in Low Gravity: A Scaling Analysis, with F. 
Rosenberger, in Low Gravity Fluid Dynamics and Transport Phenomena . 
R. L. Sani and J. N. Koster (eds.), (AIAA, Washington D.C., 1990) pp. 
87-117. (See Appendix 3) 


4. A Finite Difference Method for a Model Float Zone Free Surface 


Problem, with Y. Q. Zhang, to appear in the Journal of Numerical 
Methods in Fluids, 1990. (See Appendix 4) 
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5. The Sensitivity of a Non-lsothermal Liquid Bridge to Residual 
Acceleration , with Y. Q. Zhang, to appear in the Proceedings of the 
International Union of Theoretical and Applied Mechanics Symposium on 
Microgravity Fluid Mechanics, Bremen, September, 1991. (See Appendix 5) 

6. The Sensitivity of a Non-lsothermal Liquid Bridge to Residual 
Acceleration, to appear as an extended abstract of the above article in 
Microgravity Science and Technology, August/September 1991. 

A paper describing our work on 2D thermocapillary flow is in preparation. Results of 
work for the VIT project have been (will be) presented at the following conferences and 
workshops and invited lectures: 

1. Sensitivity of Liquid Bridge and Floating Zone Shapes Subject to Axial 
Vibrations, presented at the Alabama Materials Research Conference, University 
of Alabama in Huntsville, Sept. 20th and 21st 1989. 

2. Sensitivity of Thermocapillary Flow to Residual Acceleration, presented at 
the Alabama Materials Research Conference, University of Alabama in 
Huntsville, Sept. 20th and 21st 1989. 

3. Experiment Sensitivity: Determination of Requirements for Vibration 
Isolation, Proceedings of the NASA Workshop on Vibration Isolation 
Technology Microgravity Science Experiments, Sept. 28-29 1988. 

4. Sensitivity of Crystal Growth Experiments to Residual Accelerations, 
presented at the Gordon Conference on Gravitational Effects in Materials 
and Processes, July 30 -August 4, 1989, Plymouth State College, 
Plymouth, New Hampshire. 

5. The Sensitivity of a Non-lsothermal Liquid Bridge to Residual 
Acceleration, with Y. Q. Zhang, to be presented at the IUTAM 
Symposium on Microgravity Fluid Mechanics, Bremen, September, 1991. 

3. Liquid Bridge/Floating-Zone Sensitivity 

This section refers to Task 4 and concerns the sensitivity of an isothermal liquid 
bridge to axial acceleration and was been extended to non-isothermal situations covered 
under Task 3. The first annual report contained a detailed description of a one-dimensional 
model and the sensitivity results obtained using that model, this work was subsequently 
been refined published in the journal Physics of Fluids. A reprint is attached in Appendix 
2 . 

The 1-D model had been used to examine the sensitivity of an isothermal liquid 
bridge to axial vibration. The results indicate that the zone is most sensitive to accelerations 
with frequencies close to or equal to the lowest natural frequency of the zone. For the 
purposes of this project it is most useful to assess the sensitivity of the zone in terms of 
predicted Space Station and/or Spacelab environments. For the cases examined, we have 
seen that frequencies around the 0. 1 Hz range appear to be the most sensitive. The low 
frequency (< 10' 2 Hz) acceleration environment predicted for the space station 1 should not 
exceed levels of 10' 5 g. Higher frequencies can be associated with acceleration magnitudes 
of up to 10' 2 g. In terms of these predicted levels, or those measured on past missions 1 , 
the practical sensitivity range is restricted to disturbances with frequencies ranging from 
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10- 1 - 10 Hz. 

During the second year, work on a full axisymmetric calculation started. A model 
was developed for both isothermal and non-isothermal bridges. The adaptation of two 
numerical algorithms for the solutions of these problems has been successfully completed 
and preliminary calculations have been made. 

The early stages of development and verification of the isothermal code (and the 
ability of the method to handle buoyancy-driven flows in fixed geometries) was outlined in 
the third semi-annual report and an updated account is given below. 

For the solution of the axisymmetric isothermal and non-isothermal problems our 
objectives were: 

1. To develop a full numerical model of the axisymmetric free surface 
response of an isothermal liquid bridge (suspended between rigid disks of 
equal radius) subject to axial vibration and to compare it to the sensitivities 
predicted by a previously developed 1-D model 

2. To develop a full numerical model of the axisymmetric free surface 
response of a non-isothermal liquid bridge subject to axial vibration. The 
presence of a temperature gradient gives rise to a thermocapillary flow. We 
wish to investigate the interaction between the free surface motion induced 
by vibration, and the surface-driven flow. Interaction of these flows with 
internal buoyancy-driven flow driven by the vibration is expected to be 
negligible but is included in the basic model. 

Any algorithm that is capable of solving the equations governing the second 
objective must be capable of solving the first objective. Added complications result from 
the fact that the coupled heat transfer and Navier-Stokes-Boussinesq equations must also be 
solved to satisfy the second objective. In addition, it is convenient to use steady-state 
solutions for thermocapillary convection as initial conditions. As a result we have 
developed two algorithms which are based on the same numerical methods; one is used to 
compute steady-state velocity and temperature fields and free surface shape for the 
following problem, the other is used to calculate the time dependent response. 

3.1 Steady State interface shapes and thermocapillary convection 

Consider a cylindrical liquid zone held between two parallel coaxial circular rigid 
disks (radius = R 0 ) separated by a distance L. The liquid is a non-isothermal 
incompressible Newtonian fluid. The bridge is held between the disks by surface tension. 
The free surface of the bridge is a gas-liquid interface and is described by r=R(z,t). The 
two disks are maintained at a constant temperature T 0 . Surface heating is provided through 
a parabolic ambient temperature Too(z) and the heat transfer coefficient at die free surface is 
denoted by h. In addition, we make the assumptions that the gravitational acceleration is 
parallel to the cylinder axis and that the velocity, temperature and free surface deformation 
are axisymmetric. Furthermore, we let the surface tension at the free surface vary linearly 
with temperature and assume that the Boussinesq approximation holds. 

The governing equations are made dimensionless by scaling length, time and 
velocity with Rq, Ro/U* and U*> respectively . Here U* is a characteristic velocity given by 

IyIat 

n 

where AT=T ma x-T m i n represents the maximum temperature difference at the surface, lyl is 
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the absolute value of the derivative of the surface tension with respect to temperature, and (i 
is the dynamic viscosity. We shall refer to a "half zone" model when the ambient 
temperature has extrema at the end disks and a "full zone" when the temperature maximum 
occurs between the disks. For a full zone we shall take T max to be Too(0), and T m i n to be 

T 0 , where A = L/Ro is the aspect ratio. 

The non-dimensional pressure is 


p* + pogz 

poU * 2 


Ro 


where p* is the dimensional pressure, g is the gravitational acceleration, z is the 

dimensionless axial coordinate, and po is the density corresponding to the reference 
temperature. The temperature is rendered dimensionless using T max - T m i n . With these 
scales the dimensionless equations in a cylindrical coordinate system can be written as 

r^r + lr =0 ' (1) 


„ + w = . 3 r + _L (£“ + 1 §“ + * 2 . . jl' 

dr dz dr Re \d r 2 r dr dz 2 r 2 j 


( 2 ) 


du 

u — + w 
dr 


dw 

dz 


dp i /d 2 w , 1 dw | d 2 w\ | Gr T 

dz Re \ r dz 2 / R e 2 


( 3 ) 


dT dT i (d 2 T i dT dT \ 

dr dz Ma^dr 2 r & dz 2 ) ’ 


( 4 ) 


where the Reynolds number, Re, Marangoni number, Ma and Grashof number Gr, are 
respectively 


, M a = b^, Gr = I^. 

V (J.K v 2 

Here, v is the kinematic viscosity, k is the thermal diffusivity, (3 is the volume thermal 
expansion coefficient, and g is the gravitational acceleration. At the disks, the boundary 
conditions are 

u = w = T = 0, at z = ± y. (5) 


The symmetry conditions at the centerline r = 0 are 

dw dT 
U dr dr 


( 6 ) 
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The free surface is located at r = R(z) where the boundary conditions are: 



(7) 



jdu 

|3z 


3w\ „3R/3u 3w\ [ /3R\ 2 ]2/3T 3R3T1 

3r) + aziar'azr'i lazjj laTaz arj* 


3R _ 

u + w = 0 

a z 


( 8 ) 

(9) 


1 


'*(! f 


3T 3R 3T\ r. /'T T n 

ar B,(T - T_) = °, 


(10) 


where 


C 0 = ^l, Bi = ^, G=^| 

Yo K u 2 

and are, respectively, the capillary number, Biot number, and dimensionless gravitational 

acceleration, y 0 is the mean surface tension, and k is thermal conductivity. The force 
balance at the free surface in the normal and tangential directions are given by eqs. (7) and 
(8), respectively. Equation (9) is the kinematic boundary condition at the liquid-gas 
interface. The thermal boundary condition at the interface is given by equation (10) in 
which the equivalent heat transfer coefficient h contains the effect of the radiant and 

convective heat transfer. The constant X in (7) represents a dimensionless reference 
pressure difference 2 * 3 across the interface. In model float zone systems with fixed rigid 

endwalls, such as the one discussed here, X is determined by the following constant 
volume constraint 


V = 



K 


R 2 (z)dz = V 0 = constant . 


( 11 ) 


Finally, the condition that the contact lines between the liquid end disks are fixed is 
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R = 1 at z = ± A/2. (12) 

Numerical Method 

In the case of a two dimensional axisymmetric flow, the governing equations can be 
simplified by introducing the stream function, \| i, and vorticity, co, as new dependent 
variables: 


u -l®5t v _i 
r 3z ’ rjr 

(13) 

3u dw 

(14) 

03 ~ Ih'lh ' 


From (13), (12), (2) and (3) the following equation is obtained for co 


dco dco 

u aF + w aI 


cou _ 1 f d 2(l) 

r Re (dr 2 


, 1 ^co yco\ 1 co Gr dT 
r dr dz 2 1 Re r 2 Re 2 dr 


Substitution of (13) into (14) yields 


(15) 


rco = 


d> 

dr 2 


1 dy d V 
r dr dz 2 


(16) 


The original set of three equations governing mass and momentum has been reduced to 
two equations governing the steam function and vorticity. 

The steady free boundary problem for cylindrical liquid zone is solved iteratively, 
since the location of the free surface is a priori unknown. To obtain a solution we adopt a 
Picard iterative procedure 4 as follows: 


1. Guess the free surface shape for the initial iterate; 

2. Obtain the approximate temperature and velocity fields by transforming the governing 
equations and boundary conditions to a circular cylindrical domain via a non-orthogonal 
transformation and solve them using a pseudo-unsteady semi-implicit method; 

3. Obtain the pressure at the free surface by integrating the transformed momentum 
equation; 

4. Use the normal force balance condition at the free surface to decide how to update the 
free surface location; 

5. Return to step 2. Repeat until convergence is obtained. 

Non-orthogonal Transformation 

The region occupied by the liquid zone is transformed into a fixed circular 
cylindrical computational region using a non-orthogonal coordinate transformation, i.e. 


q = z 



(17) 


It then follows that 
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a_ j_a_ _a__c_ 3 Ra_ 

5r R ’ dz aq R a£ 


(18) 


The advantage of this transformation is that the free boundary coincides exactly with a 
coordinate line in a computational grid and regeneration of mesh during the outer iteration 
is avoided. For the sake of simplicity, a pseudo-unsteady method in association with a 
semi-implicit discretization with respect to time is used. 

The resulting system of governing equations is solved as the time (pseudo) 

derivatives of T, \\r and co — > 0. The discretization in time employs an explicit Adams- 
Bashforth scheme for the nonlinear convective terms and the implicit scheme for the 
viscous terms. Other terms in the equation are treated explicitly. For spatial discretization, 
central differences are used. The A.D.I. form for the vorticity equation is then 


where 


- 2 - (cD n+1/2 - co n .) + 2 - (L 2 co) n . - 1 (L2G))! 1 : 1 - A 
Ax \ y w 2 ' ; 'J 2 >J Re 


w 


Wk 


j. 

Re 


1 a 2 0) 5co\ n+1 , , n 

A — r + B 1 + (S^ = 0 , 


\ dn ac/ij 

- 2 - (co n . +1 - o n+1/2 ) + 2 - (L 2 tof . - 1 (Uco )": 1 - ± f 
At l u y ' 2 V ; >j 2 V y 'J Re\dn 2 /ij 


Re 


/ -.2 -> \n+i/2 

A^ + B^ +(S,r = 0 
V a; adi 


L 1 f?v a ay a \ 

1 r 2 £ \3n a£ acanj’ 


l 2 = 


’^1 AAi _u_ 
r 2 C (an ac a; an] r; ’ 


s, =- Ac 


3 2 T 1 

Ma “ aqac ’ 2 Re 


Cd c a 2 co \ | Gr 1 dT 
R 2 C dridy Re 2 R aC ’ 


09) 


( 20 ) 


S 3 = 



aqac 


+ R£w . 


Here the superscript n+1/2 denotes the intermediate step associated with the ADI 
method 7 . The velocities (at the nth step) are taken from the values at the previous step 
values. Thus, except for the thermal condition at free boundary, these are Dirichlet 
boundary conditions. The heat transfer equation is also discretized in an ADI form. The 
resulting system of discretized equations is solved using a factorization method 5 . 
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We define a steady state to occur reached when the residuals (dco/dx, 3T/3x, 3y/dx) 
of the vorticity, energy and stream function equations are less than 10' 7 . That is. 


max 



Ax 


< 10- 7 . 


(21) 


Where F represents the vorticity, temperature and stream function at the n+lth and nth 
iterative step. The numerical solution is second order accurate in space. 

Having computed the vorticity and stream functions for a given surface shape, it 
remains to iterate on the condition for the balance of force normal to the surface in order to 
obtain the final steady surface shape. In addition, the shape must satisfy the volume 
constraint (1 1) and boundary conditions. 

Two iterative schemes for determining the interface are discussed in this section. In 
both schemes, the new velocities and temperatures are taken from the current calculated 

values and the pressure at interface can be obtained by taking the T] component of 

momentum equation and integrating with respect to rj. This yields 


P = 


i fdy d w dy 9w \ i 

R \ ds dt; ds ) Re 


'a 2 w +A ^w +B 3w + c a 2 » 


, + -Or_ T 
3s3C I Re 2 ) 


jaRfr . 

Mar)) ds - 


as 2 3£ 2 a; 

The first scheme (I) involves successive approximation by the direct solution of 


, 2 Re- 1 l 3u Id R\ 2 /3w j 0R 3w \ 3R / 1 3w 3u \ 9R du 

P " 11 + "i /§Rf | R ^ ”^\ R ^ ^n’ R ^ia^jJ 

u 


= Re' 1 (q, 1 - t) r 1 




3RP 

an 


m 


ffl 2 

+ (an) a 2 R 


R 


8n 2 


(22) 


This scheme is based on the following principle. A shape is assigned to the free surface with the 

calculated pressure, velocity and temperature. An initial guess for the pressure constant X is made. 
The new interface shape is then determined directly from (22). The integral (1 1) is then evaluated 
to check whether the volume constraint has been satisfied. If it is not satisfied an inner iteration is 

made using a Newton-Raphson procedure to calculate the following improved estimate of X: 


X 


k+l 




(23) 


where 
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M/2 

AV = 7iR 2 dTi - V 0 . (24) 

J- A/2 

The above procedure for determining X is quite effective and is repeated until the volume 
constraint is satisfied. RCn) is then updated. New velocity, pressure and temperature fields 
are calculated using the updated value of R(tj). The outer iteration is repeated until 

max |R n+1 - < e , (25) 

where we take e = 10* 4 . 

The second scheme (II), used by Ryskin and Leal 7 , uses the residual of the normal 
component of force balance to drive the shape to its steady position. This is equivalent to 
equating the residual with an artificial capillary force. This effective force causes a local 
displacement of the interface in the direction of the force. The magnitude of the local 
displacement is proportional to this force. The interface shape at each iteration is thus 
modified so as to reduce the residual until the condition (22) is met. It follows that at each 
iteration the improved interface shape is given by 

Rj n+1 = R m + aExj, (26) 

where Exj is the residual of the force balance equation at the jth surface location and the 
constant coefficient a is determined by numerical experiment. In order to ensure 
convergence, a should be small. We found that the values of a which led to rapid 

convergence depended on the product of Re' 1 and Co' 1 . If oc is chosen to be too small the 
amount of CPU time used increases subst antial ly. 

The change in volume between the mth and (m+l)th iteration can be found from the 
volume constraint (11) and equation (26) and neglecting higher order terms, i.e. 

M/2 

Exj R™ dr\ = 0. (27) 

J A/2 

The pressure constant X is contained in Exj and is obtained by satisfying (24). Even then 
the liquid bridge may still change volume slightly at each iteration owing to numerical error 
and higher-order effects. These small changes can accumulate and eventually result in a 
gross error. To prevent this, the formula (26) is modified to 

R m + i = R m(V SL J 1/2 + aEx . . (28) 

\V m / 

We have used the method described earlier to examine the influence of various 
parameters on momentum and heat transport and meniscus shape. In addition, a 
comparison between the results obtained with the combined scheme and those of Hyer et 
al. 7 has been made. 

A parametric study of the effects of varying the temperature difference AT (keeping 
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Ma/Gr fixed), Marangoni number, Ma, Nusselt number, Nu, aspect ratio. A, and Grashof 
number, Gr, has been made and are described in Appendix 4. 

Unsteady interface shapes (g-jitter response) for isothermal and non-isothermal 
bridges. 

To obtain the time dependent internal and free surface response of isothermal and 
non-isothermal liquid bridges to axial acceleration we employed the following procedure. 

To obtain a solution we follow the approach of Kang and Leal 8 and adopt a Picard iterative 
procedure 4 as follows: 

1. Use the free surface shape, stream function, vorticity, velocities, and temperature 
computed from the steady calculation as the initial condition. (For the isothermal case the 
liquid bridge is initially static with a free surface shape given by the prevailing steady axial 
acceleration ; 

2. To obtain the temperature and velocity fields at the next time step we transform the 
governing equations and boundary conditions to a circular cylindrical domain via a non- 
orthogonal transformation and solve them using a semi-implicit Adams-Bashforth/Crank- 
Nicolson time discretization. Centered differences are used for the spatial approximation ; 

3. Obtain the pressure at the free surface by integrating the transformed momentum 
equation; 

4. Use the normal force balance condition at free surface to determine the updated free 
surface location at the new time step; 

5. Return to step 2. Repeat until convergence is obtained. 

Our early results showed good agreement between the predictions of the 1-D 
isothermal model and the axisymmetric calculations (see Table 1) Further calculations have 
been made and for the cases examined to date these findings are confirmed. For the cases 
examined to date, the effect of thermo-capillary convection on the response of the free 
surface shape has been found to be negligible Our results are given in Appendix 5 

Table 1 : Comparison of Bridge Radii (R) (to 3 decimal places) at two locations predicted by 

1-D and axisymmetric (2D) models, g= 1.42x1 O' 3 ms -2 , the frequency of axial vibration is 
3 Hz. 


Time (s) 

0 

0.22 

0.39 

0.48 

Ri-d 

l.OOl 

1.149 

1.017 

0.891 


1.001 

1.146 

0.998 

0.919 

|Ri-d 

1.001 

1.159 

0.994 

6.909 


1.001 

1.150 

0.986 

0.929 


4. Thermocapillary Convection -2D 

The work for Task 2 involves the solution of the incompressible Navier-Stokes - 
Boussinesq equations in a rectangular region with one "free surface" boundary condition 
which equates the tangential viscous force at the surface to the temperature gradient along 
the surface. A pseudo-spectral Chebyshev collocation method 5 has been employed. The 
basic approach and the equations to be solved were outlined in the first annual report. 
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The objective is to study the effect of g-jitter on thermocapillary flows for a range of 
dynamic Bond numbers Bd = pgL 2 /[dy/dT] (here the denominator represents the change in 
surface tension with respect to temperature). To date we have examined fluids with Prandtl 
numbers in the range 10' 2 to 10, which spans the range associated with most experimental 
materials. Marangoni numbers in the range 1 to 10 3 have been subjected to single 
frequency excitations with amplitudes in the range HH-IO* 2 g. As expected, the sensitivity 
of the flow to a given residual acceleration depends on the relative magnitudes of the 
Reynolds number and the Grashof number (see discussion in the 1st annual report and 3rd 
semi-annual report). 

Description of the model system 

The physical system currendy under examination involves a 2-D rectangular region 
which is bounded by three rigid walls and a "free surface". Two of the rigid walls are 
isothermal with temperatures Th and Tc- The other boundaries are adiabatic (no 
perpendicular heat flux) and the "free surface" is constrained to be a straight line. This 
retains the essential character of a free liquid surface while removing the numerical 
problems associated with a deformable one. The forces at the interface may be decomposed 
into components normal and tangent to it. Contributions to the force balance include the 
fluid pressure, the viscous forces, and the capillary forces. For the normal component of 
force balance, the difference between the pressure and viscous forces on either side of the 

interface are proportional to the interface curvature multiplied by the surface tension y. For 
the tangential force balance at the interface, the difference in the shear stress components is 
balanced by a gradient in surface tension. Such gradients will arise when there is a 
temperature gradient along the interface. The gradient in temperature produces a gradient in 
surface tension that induces shear stresses and, thus, fluid (thermocapillary) motion. In our 
calculations the residual gravity vector is capable of taking any orientation with respect to 
the boundaries and can be time-dependent. 

We approached this problem using a pseudo- spectral technique which incorporated 
the influence matrix method. Comparison of our results for steady flows with results 
obtained by Carpenter and Homsy [9] and Zebib et al. [10] indicate that our method has 
superior spatial accuracy over finite difference and finite volume methods. (See 
comparisons in tables 2 and 3. ) Our g-jitter results confirm our earlier predictions that, for 
Ma»l, if the free surface is not allowed to deform thermocapillary effects dominate for the 
range of container dimensions, acceleration magnitudes and frequencies typically associated 
with spacecraft experiments. The results of our work comprise the Master’s thesis 
(Universitie D’Aix-Marseille, France) of Ms. Helene Cordier who worked at the CMMR 
this year. A paper intended for journal publication which describes our results is now being 
prepared. 

5. Thermo-capillary Convection - Three Dimensional and Axisymmetric 

Owing to the premature termination of the grant this work could not be completed 
although some progress was made during the no-cost extension. This work is described in 
Appendix 3. 

6. Transport in Closed Containers 

In the original proposal we agreed to make available the results of related work at 
the CMMR not funded through this grant. These results were included in the appendices of 
the 3rd semi-annual report and a summary is given in the paper in Appendix 1. 
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Appendix 1 





Low-Gravity Experiment Sensitivity to 
Residual Acceleration: A Review 


In this article, work related to the analysis of experiment sensi- 
tivity to residual acceleration experienced in low earth orbit 
spacecraft is reviewed. Most of the work discussed concerns 
heat, mass and momentum transfer in fluid physics and materi- 
als science experiments . On the basis of our current under- 
standing of experiment sensitivity it is concluded that, in gener- 
al, experimenters should be concerned about the effect of 
residual acceleration and that careful modelling included as 
part of an experiment program will enable optimal use of the 
limited experiment time available in space . 


1 Introduction 

It has been recognized for some time that the low-gravity ac- 
celeration environment associated with a spacecraft in low 
earth orbit offers an opportunity to study certain physical 
processes which are difficult to investigate under the gravita- 
tional acceleration experienced at the earth s surface. Our ex- 
perience with such opportunities has led us to realize that the 
residual acceleration environment on board an orbiting 
spacecraft is not as low or as steady as would be desired lor 
certain classes of experiments. The sources of the residual ac- 
celeration include [1-5] crew motions, mechanical vibrations 
(pumps, motors, excitations of natural frequencies ot space- 
craft structures), spacecraft maneuvers and basic attitude mo- 
tion, atmospheric drag and the earth’s gravity gradient. The 
accelerations are characterized by temporal variations in 
both magnitude and orientation. Such disturbances will, in 
particular, affect those experiments susceptible to buoyancy 
induced fluid motion [6]. Indeed, recent studies [7, 9] indicate 
that transient disturbances can have undesirable long term 
effects. The analysis of the sensitivity of any proposed flight 
experiment is necessary for a variety of reasons. These in- 
clude the need to optimize the use oi the limited time avail- 
able for flight experiments, the interpretation of experimental 
results and the determination of tolerable acceleration levels 
to be used in planning for NASA's Space Station [3]. 

Examples of experiments that are conducted in the low- 
gravity environment include critical point studies [10]. crystal 
growth [11-18], and diffusion experiments [19, 20]. Since 
these experiments involve density gradients in the fluid 
phase, they are inherently sensitive to the effects of buoyan- 
cy-driven fluid motion [21-25]. For certain types of systems it 
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is known [7, 8, 26-37] that the associated physical processes 
(mass transfer, heat transfer, convection, oscillation and dis- 
tortion of free surfaces) are sensitive to time-dependent ac- 
celerations. For example, the transfer of mass and heat in 
fluid systems can be significantly affected by oscillatory 
flows [28-37]. The effect of oscillation enhanced transport 
conditions on the local variation in composition during crys- 
tal growth is not well characterized for most of the systems 
relevant to crystal growth under microgravity conditions. Or- 
der of magnitude analyses [7, 8, 38] suggest that for certain 
combinations of physical properties and growth conditions, 
oscillations in the residual acceleration may adversely affect 
the mass transport conditions. For instance, the DMOS ex- 
periment [39] on STS 51-D showed extensive evidence of 
convective mixing of liquids. The degree of mixing is greater 
than can be attributed to the quasi-steady low-gravity com- 
ponent, but can be accounted for when oscillations in effec- 
tive gravity are included in the description of the transport 
conditions. 

It is possible that in future space laboratories some experi- 
ments will require isolation from vibration. The design of ef- 
ficient isolation systems, however, requires determination of 
the tolerable acceleration levels (amplitudes and frequencies) 
for given experiments. The results of sensitivity analyses can 
be used in the development of vibration isolation systems to 
identify those experiments that will need isolation and by 
supplying data concerning frequency dependence, the effect 
of transients etc. The system can then be designed to filter 
out those bandwidths to which the experiment is predicted to 
respond adversely. 

Low-gravity experiments which involve liquids with free 
surfaces are most likely to be sensitive to vibrations. Experi- 
ments of this type include studies of liquid bridges, equilibri- 
um shapes of drops and bubbles [40-43], and thermocapil- 
lary flow experiments [43-46]. Evidence for the oscillation of 
liquid zone shapes was found on the D-l mission. Long liq- 
uid bridges were more sensitive to residual gravity, and ex- 
hibited random oscillations in zone shape [40]. 

It should be recognized that the sensitivity of a given pro- 
cess to the overall low-gravity environment will depend on 
one or more of the following factors: 

a) magnitude and direction of accelerations 

b) system geometry 

cj boundary conditions (e.g. insulated or conducting walls, 

solute sources and sinks), and 

d) physical properties of the participating materials (viscosi- 
ty, thermal and solute diffusivities). 

Any analysis of the sensitivity to oscillatory accelerations and 
transient disturbances characteristic of the anticipated envi- 
ronment aboard Space Station, or any other spacecraft, must 
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take these factors into account. In general, a full mathemati- 
cal characterization may not be practical. However, identifi- 
cation of sensitive experiments in combination with a judi- 
cious choice of simplified numerical models can be used to 
choose operating parameters that will optimize the use of the 
low-gravity environment. With this in mind, we review work 
related to the effects of steady and time-dependent residual 
accelerations, and associated work concerning transport in 
oscillating flows. 

In sec. 2 we describe the components of the low-gravity 
environment. In sec. 3 previous work involving order of 
magnitude (OM) estimates is discussed. The application of 
OM analyses to acceleration-sensitivity determination is dis- 
cussed in sec. 3.1. The results of simple analytical models 
are summarized in sec, 3.2. Work related to linear stability in 
the presence of modulated gravity is briefly discussed in sec. 
3.3. The results of direct numerical simulations are reviewed 
in sec. 3.4. In sec. 4, the results and utility of residual accel- 
eration analyses are summarized and discussed. 


2 The Microgravity Environment 

The low-gravity environment of a spacecraft has been re- 
ferred to as a "zero gravity” environment. This is because any 
object within the craft is subject to roughly the same acceler- 
ation due to the earth’s gravitational field. Thus, a tree object 
may move along approximately the same orbital path as the 
spacecraft. As a consequence, to an observer in the space- 
craft objects appear to behave as if no gravity were present. 
However, any body capable of motion relative to the space- 
craft will experience an acceleration relative to the mass cen- 
ter of the spacecraft. Tins relative acceleration arises from 
several sources which include the gravity gradient of the 
earth, atmospheric drag and attitude motions, as well as ma- 
chinery vibrations and crew activities. 

The relative acceleration that is associated with the earth's 
gravity gradient is defined as follows. The mass center of the 
spacecraft is subject to a force F(r 0 ) due to the gravitational 
attraction of the earth, where r t j is the position ot the mass 
center of the craft with respect to the mass center ot the 
earth. A panicle at a position r within the spacecraft is sub- 
ject to a force F{r). If, as a lirst approximation, we ignore the 
interaction between the spacecraft and the mass of the pani- 
cle, and if the distance i r — r§ \ is small compared to the se- 
mi-major axis of the spacecraft orbit, the difference between 
the forces applied at rand r 0 define the gravity gradient G as 
follows [2-6] 

F(r) - F(fo) = Gz, C = r- r, . , 

If the spacecraft maintains a fixed orientation, for example 
with respect to the sun. there is no rotation of the spacecraft 
frame relative to the geocentric frame. In this type ot orbit 
both the magnitude and the orientation of the residual accel- 
eration vector change as a function of time and there will be 
no steady residual acceleration. 

In addition to the gravity gradient acceleration the atmo- 
spheric drag force on the spacecraft can be significant [2-6], 
despite the fact that the atmospheric density at typical shuttle 
altitudes is only on the order of 10“ u kg m \ Atmospheric 
drag causes a slow inward spiral of the spacecratt. To an ob- 
server in the spacecraft the effect ol atmospheric drag is to 
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Table 1 Changes in residual acceleration components [ug nr 7 along 
the local vertical iXu and perpendicular to the oribital plane (X : ) asso- 
ciated with the gravity gradient stabilized attitude , and estimated ah 
moshpenc drag accelerations, a and b are, respectively, best and worst 
drag estimates for the Shuttle (I) and NASA \s Space Station (2). From 


altitude 

X\ 

X- 

drag 1 a 

drag lb 

drag 2 a 

drag 2 b 

[km] 







275 

0.414 

-0.137 

2.10 

0.360 

10.6 

3.1 

300 

0.410 

-0.136 

1.20 

0.210 

6.1 

1.8 

325 

0.405 

-0.134 

0.69 

0.120 

3.6 

1.1 

350 

0.40 T 

-0.133 

0.41 

0.070 

2.1 

0.6 

400 

0.392 

-0.130 

0.24 

0.025 

0.7 

0.2 

450 

0.383 

-0.127 

0.05 • 

0.009 

0.4 

0.1 

500 

0.375 

-0.124 

0.02 

0.001 

0.8 

0.02 


produce a relative acceleration that is equal in magnitude but 
opposite in sign to the atmospheric drag force. Table 1 gives 
atmospheric drag acceleration estimates and relative acceler- 
ation related to the gravity gradient and centrifugal force in 
ug per m from the mass center of a spacecraft that is in a cir- 
cular orbit with a gravity gradient stabilized attitude. X x is 
parallel to the local vertical and X 2 is perpendicular to the 
orbital plane. For this attitude the acceleration points away 
from the mass center along X\ and toward the mass center 

along X 2 - ' . 

If the spacecraft does not maintain such a fixed orienta- 
tion then, in addition to (1), centrifugal, Coriolis and Euler 
accelerations will become apparent in the spacecraft refer- 
ence frame. They arise in conjuction with spacecraft attitude 
motions and depend on the nature of the rotation of the 
spacecraft frame of reference relative to a Fixed geocentric 
frame. If a spacecraft in a quasi-circular orbit continuously 
rotates relative to the fixed geocentric frame such that a giv- 
en direction in the spacecraft frame is always oriented paral- 
lel to the craft’s position vector r 0 , a centrifugal acceleration 
will augment the gravity gradient in the direction parallel to 
the position vector, and cancel the component tangent to the 
flight path. Under these conditions there will be a steady 
component to the residual acceleration vector. In practice 
this situation arises, for example, in the so-called gravity gra- 
dient stabilized attitude. 

Euler accelerations arise when the rate of rotation ot the 
spacecraft frame changes as a function of time, as it would in 
a non-circular orbit, or as the semi-major axis of the orbit 
changes as the craft slowly spirals inward due to atmospheric 
drag Coriolis accelerations will also be apparent but are pro- 
portional to the velocity of an object relative to the space- 
craft. Thus, unless relative velocities are large, the Coriolis ef- 
fect will be insignificant [6]. 

In addition to the above, there are a variety of residual ac- 
celeration, components that occur over a broad range of am- 
plitudes and frequencies. The nature and sources of these ac- 
celerations have been discussed in [1-3]. Figs. 1 and 2 are 
examples of residual accelerations measured on orbit. 

3 Analyses of Residual Acceleration Effects 

As discussed in the introduction, many experiments are un- 
dertaken in orbiting spacecratt in order to reduce or elimi- 
nate unwanted effects of buoyancy driven motions in fluids. 
Sources of buoyancy in fluids include discrete interfaces, as 
they occur, for example between immiscible fluids, solid par- 
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Time {sec] 




Frequency [Hz] 



Fig. 1. An example of a relatively "quiet period " measured on SL-3 
(< 4- 10~- g. g = 9.8 m s~’) 

a) Total acceleration array, b) Combined amplitude spectrum, note 
components in the 5-6 Hz and 7-8 Hz, 10-12 Hz, 17 Hz and 30-35 
Hz ranges . A fter Rogers and Alexander [47] 


Fig. 2. “Window" containing three thruster pulses 
a) Total acceleration array, maximum magnitude ~ 1- 10~ 2 g, b) Com- 
bined amplitude spectrum, dominant frequency components 11 Hz and 
17 Hz. After Rogers and Alexander [47] 


tides and fluids, and bubbles and fluids with different densi- 
ties, other buoyancy sources are internal density gradients 
that arise as a consequence of gradients in temperature and 
concentration. The nature of the response of a bulk fluid to 
the effective body force associated with a low-gravity envi- 
ronment will depend, in part, on the relative orientation of 
the residual acceleration vector and the density gradient. In 
particular, if the integral of the effective body force Qb 
(whefe gls the density and ^ the effective acceleration) acting 
on any closed circuit within the fluid is non-zero, motion will 
immediately ensue. That is, if 

^ d c ■ Qb — ^ d4 • curl Qb 

c s 

= \ dA * fgrad^x b + pcurl 6] # 0 , (2) 

motion will occur. Note also, that even if there is no density 
gradient motion will occur if 


\ dA ■ curl Qb ± 0 . (3) 

s 

This includes the case of “precessional stirring” [48]. This 
could occur, for example, when a spacecraft orbiting in a ra- 
dial attitude (i.e. it continuously rotates about an axis per- 
pendicular to the orbital plane) underwent an additional ro- 
tation about some axis other than the basic orbital axis. 

There have been three approaches to the determination 
experiment sensitivity to the residual acceleration environ- 
ment. 

a) order of magnitude analyses [38, 49-52], 

b) analytical models [5, 6, 39, 53-59], and 

c) direct numerical simulation [9, 3, 60-73]. 

In addition to work directly related to microgravity, there is 
extensive literature on mass and heat transport in oscillatory 
flows [26-36], the effect of gravity modulation on the classi- 
cal Benard problem [74-77] and the linear stability of planar 
fluid surfaces [78, 79]. 
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3. 1 Order of Magnitude Analyses 

An r/rder of magnitude (OM) or dimensional analysis of a 
physical process involves the examination of the system of 
equations that are assumed to govern the process. For fluid 
physics and material science experiments, the relevant physi- 
cal processes are governed by equations describing the trans- 
port of heat, mass and momentum in single and multi-com- 
ponent fluid systems. The usual approach to OM analyses 
[51] involyes the definition of reference quantities which ap- 
propriately characterize the physical system (i.e. reference 
scales for velocity, time, length, forces, etc.). Application of 
these scales to the (dimensional) governing equations, fol- 
lowed by the definition of characteristic dimensionless 
groups admits a comparison of the order ot magnitude of 
every term in each equation. There are several comprehen- 
sive discussions of OM analyses applied to low-gravity situa- 
tions (e.g. [38, 49-51]). 

The advantage of the OM approach is that a great deal of 
information pertaining to the sensitivity of a given experi- 
ment can be obtained with little computational effort. The 
disadvantages are related to the fact that the choice of char- 
acteristic reference quantities for length, velocity, time etc. 
are not, in general, known a priori [80]. Thus, in cases for 
which the chosen characteristic scale is not (for a given set of 
conditions) representative of the actual process, the results of 
the analysis may be in error, often by several orders of mag- 
nitude. In addition, application of OM analyses to the prob- 
lem of g-sensitivity has invariably involved implicit lineariza- 
tion and often poorly represents the multi-dimensional 
nature of the physical process. This can also lead to incorrect 
OM predictions. Finally, analyses are generally restricted to 
the examination of a single component disturbance (e.g. 
- sin(tu/)). As a result, they may not properly indicate the 


response of the system to the typically complex multi-fre- 
quency disturbances characteristic of the spacecraft environ- 
ment. It is nonetheless useful to examine the general trends 
predicted by such analyses since they can provide at least 
qualitative guidance for more detailed analytical and numeri- 
cal studies. 

The determination of tolerable g-Ievels using OM analyses 
is based upon estimations of the response of a system to ei- 
ther steady accelerations or time-dependent disturbances of 
the form 

g(t) = acos(cot), W 

where co is the angular frequency of the disturbance. The sus- 
ceptibility of an experimental system to such disturbances is 
defined via the magnitude of a particular response (say a 
temperature, compositional or velocity fluctuation) which 
must not be exceeded in order to bring the experiment to a 
successful conclusion. The most obvious trend predicted by 
analyses to date is shown in figs. 3-5. The curves depict the 
maximum tolerable residual acceleration magnitude as a 
function of frequency ,/ = co/2 n, for given experiments and 
have the form 

8 max = FfW, 77 ( 5 ) 

where fl cnt is the magnitude of the maximum allowable re- 
sponse, and the p h i - 1, N represent the N material proper- 
ties of the system (such as viscosity, thermal and species dif- 
fusivities). In general, below 10“ 2 to 10~ 3 Hz the response is 
expected to be either independent or weakly dependent on 
frequency. At higher frequencies the analyses indicate that 
the allowable residual acceleration magnitude increases as 
approximately the square of the frequency. This behavior can 





Fig. 3. Estimated tolerable residual accelerations (in units of g — 
9.8 m s~ 2 ) as a function of frequency for: 1) a fluid physics experiment 
involving a temperature gradient, 2} a crystal growth experiment and 3) 
a thermodiffusion experiment. After Monti etal. [38] 


Fig. 4. Estimated tolerable residual accelerations (in units of g = 
9.8 m s~ 2 ) for semiconductor and metal solidification experiments. The 
sensitivity parameter is longitudinal segregation, Ac/c. Pe s =* 5 , tem- 
perature gradient — 50 K cm~ } , k = 0.1. After Monti etal. [36] 

Fig. 5. Estimated tolerable residual acceleration (in units of g — 
9.8 m s~ 2 ) ranges for diffusion experiments. After Monti etal. [36] 
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he understood upon consideration of the following OM esti- 
mates [38] (used to obtain fig. 3) for velocity, temperature and 
concentration changes associated with the acceleration (4): 


r, - 


ST- K 


SC - K 


£144 L£_ 

(co : + 2 

AT L 


( co 2 + k-'-L*)''' ’ 

AC L 

(a 2 + D : /L 4 ) 12 ’ 


( 6 ) 

(7) 

( 8 ) 


where g*is the residual acceleration magnitude, 
n is the mass density, 

and Aq/q\s the change in density due to temperature and/or 
composition; 

.A T/L, AC L and L are. respectively, a characteristic temper- 
ature gradient, concentration gradient and length; 
k is the thermal diffusivity. 

The sensitivity limits are obtained by using either ST, SC 
or V g as a sensitivity parameter and setting g* equal to the al- 
lowable residual acceleration g max . 

For situations which involve mixed thermocapillary con- 
vection and buoyancy-driven convection the governing di- 
mensionless groups are [51]: 
the Bond number, Bo — ng* Lr/y , 
the Grashof number Gr = AT ' ft T g* LVv 2 , 
the Prandtl number Pr = v/k % 

and the Surface Reynolds number Rs = |y r ! AT L/liv. 

Here ! y T \ is the rate of change of surface tension y with tem- 
perature, 

jl r is the coefficient of thermal expansion, 
ris the kinematic viscosity, 
u is the shear viscosity. 

These groups respectively represent the ratios of buoyan- 
cy force to surface force, buoyancy force to inertial force, 

* thermal and momentum timescales, and surface force to iner- 
tial force. 

The relative importance of gravity and thermocapillary 
forces can then be estimated upon considering the relative 


magnitudes of the velocity due to buoyancy and the velocity 
due to thermocapillary effects [50, 51]. The latter is given by 

v* = IzlLAI . ( 9 ) 

For thermocapillary flow to dominate, the ratio of V g to V * 
must be less than one. For this to occur, the maximum mag- 
nitude of the acceleration must satisfy: 


g* < i yt\ 


{of ± u 2 /L 4 ) i/2 
vPt 


( 10 ) 


Eq. (10) can be re-interpreted in terms of the ratio of the 
Grashof number, Gr, and the surface Reynolds number, Rs. 
The condition for buoyancy and thermocapillary forces to be 
the same magnitude is 

Gr =4ts{Sf~ 4- \) ] 2 , (II) 

where St = ojL 2 / vis the Strouhal number. 

This condition is illustrated in Ftg. 6. Three regimes of in- 
terest are identified. These regions are defined by the relative 
values of the Grashof, Strouhal and surface Reynolds num- 
bers. In region I, the Strouhal number, St, is less than one, 
and the condition that thermocapillary forces dominate is 
equivalent to that for a steady flow. In region II, St is greater 
than one. In this region, for a Fixed value of Rs t the value of 
Gr required to give buoyancy forces equal weight to thermo- 
capillary forces increases with increasing values of the Strou- 
hal number. This essentially reflects the fact that the charac- 
teristic time for the fluid response greatly exceeds the period 
of the disturbance. Thus, thermocapillary forces dominate at 
higher values than in region I. In region III buoyancy forces 
predominate over thermocapillary forces. For situations here 
a boundary layer scaling is appropriate i.e. Rs > 1 [50]. The 
relative importance of buoyancy can be estimated via the ra- 
tio 


Gr = RF 5 (S/ : 4- l) 1 : . (12) 

Fig. 7 depicts the relationship between the estimated toler- 
able gdevel and the frequency of the acceleration for experi- 



sui 


6 Strouhal Number (St) 



Fig. 6. Tolerable acceleration limits for a thermocapillary flow experi- 
ment expressed in terms of Grashof Gr l, Strouhal fSf and Surface 
Reynolds iRs) numbers 


Fig. 7. Tolerable acceleration limits fin units of g =«• 9.8 m s~ 2 ) as a 
function oj frequency for the thermo-capillary’ experiments with proper- 
ties listed in table 2. aj Gravity forces I°o of thermo-capillary' forces, 
b) gravity forces I0°6 of thermo-capillary' forces 
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*e 2. Physical properties and dimensionless groups for various experimental systems: I -di-methyl silicone oil [44], 2-methanol [81 ], 3-NaNOi [45], 4-indium 
5-selenium - [82], 6-silicon fl5] 
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ments with the physical properties given in table 2. These fig- 
ures represent estimates based on two extreme How regimes, 
penetrative or “viscous” flow (fig. 7 a) and boundary layer 
flow (fig. 7b) and may be respectively considered as “worst” 
and “best” case estimates for the given sensitivity parameters. 
The curves have been determined for a sensitivity criterion 
which requires the tolerable acceleration levels to be such 
that the buoyancy forces should not exceed 10% of the ther- 
mocapillary forces (if the requirement had been 1 %, the tol- 
erable accelerations would be an order of magnitude small- 
er). In terms of acceleration levels predicted for the space 
station or measured on past missions the practical sensitive 
ranges are restricted to disturbances with frequencies less 
than 10 Hz. The majority of the estimates suggest that for 
these experimental parameters, periodic vibrations will gen- 
erally not lead to significant buoyancy effects. 

The above examples estimate the effects of simple single 
component accelerations. In reality, low-gravity disturbances 
tend to be associated with more than one frequency. Given 
that a system may respond to a multi-frequency disturbance 
in an “additive” way, tolerance curves such as these may un- 
derestimate the response to a given low-gravity disturbance. 
For example, models of the DMOS experiment, show that 
mass transport has an additive response to the (multi-fre- 
quency) residual acceleration [39]. Post flight analysis of the 
experiment results has demonstrated that the amount of mix- 
ing observed between organic liquids can be explained by 
the additive response of the system to a multicomponent dis- 
turbance. 

In many cases an a priori choice of length or time scales 
may not be obvious. For example, in a study of dopant (so- 
lute) uniformity in directionally solidified crystals, the OM 
estimates of Rouzaudz tal. [52] and Camel and Favier[ 83] are 
in agreement with the direct numerical simulations of Chang 
and Brown [84] for a Schmidt number (ratio of the melfs 
kinematic viscosity to dopant diffusivity) of 50 but overesti- 


mate the amount of radial segregation for a Schmidt number 
of 10. 

The results of an order of magnitude or scaling analysis 
have been compared in [80] with those of numerical simula- 
tions of the effects of steady residual acceleration on compo- 
sitional non-uniformity in directionally solidified crystals. 
The basic model consisted of a Bridgman-type system with 
an ampoule translated between hot and cold zones spearated 
by a thermal barrier. In order to apply the CameFFavier 
technique [52, 83] to the numerical simulation, it was neces- 
sary to evaluate a proportionality constant, ,r 0 : . This relates 
the maximum convective velocity obtained in the numerical 
simulations with the Grashof number. The value of this con- 
stant depended on the orientation of the acceleration, and on 
the chosen length scale. The adiabatic zone length L a yielded 
results which best showed the trends predicted by the Camel- 
Favier approach. Table 3 lists values of Xq associated with 
different choices of L a as a reference length. Having calculat- 
ed .Vo we then determined the transport regime (see refs. [52, 
80, 83] for details) by graphing, in fig. 8, GrSc vs Pe g for the 
case of the residual acceleration parallel to the interface. 
Here Sc (Schmidt number) gives the ratio of kinematic vis- 
cosity to thermal diffusivity, and Pe g — VL/D , is a ratio of 
the translation rate to the characteristic diffusive speed where 
L is the characteristic diffusion length and D the solute diffu- 
sivity. The lateral composition non-uniformities in the crystal 
obtained from numerical simulations [9] are also given in 
fig. 8. A comparison between computed and predicted non- 
uniformity values showed that at high values of GrSc (region 
lc in tig. 8), the order of magnitude approach meets with 
some success but is less faithful at low values of GrSc (for ex- 
ample the transition between regions lb and II). Fig. 8 also 
contains two cases for which Pe g is reduced at fixed GrSc. 
According to the estimates the lower GrSc case should have 
yielded a higher compositional non-uniformity as the growth 
rate was reduced, i.e. the estimates predict that the system 


Table 3. Comparison of numerical (2’j and order of magnitude estimates („,) of compositional non-uniformity for different choices oj length scale 
[80 j. The non uniformity x in these cases is given by = = \ c(0) - c x ; c(0). c* is the melt composition far from the crystal and, because k < l,c(0) 

is the maximum value of the composition at the interface 
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Fig. 8. A plot of GrSc vs. Pe.for the results of numerical simulations for 
steady accelerations parallel to the crystal-melt interface [80]. The solid 
and dashed lines delineate the following transport regimes: l a-c con- 
vective regime; 11 advecti ve-d iffusi ve regime [S3]. Note that the adia- 
batic zone length has been used to calculate the Grashof and growth 
Peclet numbers, and that x t) = 0.27. Numbers in parentheses indicate 
the value of compositional non-uniformity 2 ~ c WHtx - c, ltlw fc\ ar • 
100%, wehre t\ denotes the crystal composition at the melt-crystal inter- 
face, associated with each of the numerical results 


moves closer to the convective diffusive transition. This was 
not reflected in the non-uniformity obtained from the numer- 
ical simulations. Evaluation of the numerical results showed 
that the increase in the characteristic diffusion length scale 
(associated with reduction in growth rate) led to an order of 
magnitude decrease in the solute gradient at the interface. 
Hence, because the system-wide variation in composition is 
small, the interfacial non-uniformity is correspondingly 
small. 

The above results indicate that scaling arguments alone 
can at best be expected to yield order of magnitude accuracy. 
The limitations of scaling arguments are hardly surprising if 
one considers the neglect of the multi-dimensionality 
(boundary conditions, How structure, etc.) of the physical sit- 
uation which is inherent in such scalar descriptions. Of 
course, prior knowledge of the system behavior can be used 
to locally improve the accuracy of the scaling, (e.g. Kimura 
and Bejan [85] ). Similarly, if a system locally exhibits one-di- 
mensional behavior, such as boundary layer flow, then the 
appropriate scaling for the formal reduction of the transport 
equations can be used effectively to estimate the local system 
response. In general, however, our understanding of the re- 
sponse of systems to residual acceleration can be furthered 
only by an approach which uses scaling, mathematical mo- 
delling and, naturally, the results of experiments in a comple- 
mentary fashion. 


3.2 Analytical Models 

A few attempts have been made to assess the effects of time- 
dependent residual acceleration using analytical models 
[53-59]. The first attempt to model the effects of disturbances 
on heat transfer between a fluid and a solid bounding medi- 
um was carried out by Gebhard [53]. Under the assumption 
that the fluid and its container would be subject to a “se- 
quence of abrupt relative displacements spaced by a time in- 
terval r f . . Gebhard took r c to be a random distributed 
variable. He showed that as r m (the most likely value of r f ) 
decreases, the resulting heat transfer was greater than pre- 
dicted by conduction in the absence of a disturbance. In par- 
ticular, it was found that for r m < d 2 /K (where d is a charac- 
teristic length and k is the thermal diffusivity) the heat 
transfer increased rapidly as r m decreased. The main limita- 
tion of this model is that fluid convection (and thus fluid 
properties such as kinematic viscosity) do not enter the ana- 
lysis. It is well known that the response of heat transfer to 
convection varies according to the Prandtl number. In partic- 
ular, for low-gravity flows (as we shall see later) the tempera- 
ture field is less sensitive to convection for low Prandtl num- 
ber fluids ( Pr < 1). 

Approximate solutions for transient convection in a cylin- 
der with an azimuthal temperature variation have been ob- 
tained by Dressier [54]. Residual accelerations representing 
the motion of an astronaut and a transient rotation of the 
spacecraft were imposed on the system. Linear accelerations 
were imposed perpendicular to the cylinder axis. For kine- 
matic viscosities of 10~ 2 cm 2 s _I , cylinders of radius 1 and 
2 cm with a maximum temperature difference of 95 K across 
the diameter were examined. The maximum fluid velocities 
consequent to the astronaut motion (modelled by a sequence 
of two step functions separated by 2 s with 4 10 4 g magni- 



Fig. 9. Transient convection in large and small test cells due to residual 
accelerations caused by astronaut motion after Dressier [54]. Here x is 
the angular displacement, v the velocity and t ; j is the decay time given 
by 0.68 ■ R : v 
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-tude, 0.5 s duration and zero mean, oriented perpendicular to 
the temperature gradient) were 12 pm s _l for both the large 
and small diameter cylinders (see fig. 9). For the larger cylin- 
der the second pulse reduced the velocity from 12 um s -1 to 
0 by the end of the pulse. The effect of the second pulse was 
to change the velocity in the small container from 12 to —3 
um s' 1 by the end of the pulse. The velocity decayed asymp- 
totically to zero reaching a value close to zero after 15 s. The 
reason the response was different for the two containers can 
be explained simply in terms of the momentum decay times 
(t ],2 ~ R 2 /v) for each container. The velocity in the larger 
container had hardly begun to decay when the second im- 
pulse arrived. The shorter decay time associated with the 
smaller container meant that the impulse acted on a smaller 
magnitude velocity which explains the change in direction of 
the motion. While the disturbance had a mean value of zero 
the fluid particles attained a finite displacement (indicated by 
the value .v of the angular displacement in fig. 6). Similarly, a 
disturbance corresponding to a simple spacecraft maneuver 
was examined. The simulated maneuver corresponded to a 
90° rotation of the spacecraft at a rate of 1° s _l . The maxi- 
mum velocities associated with the centrifugal acceleration 
were, respectively, 40 and 35 pm s~’ for the large and small 
diameter cylinders. 

Experiments involving free surfaces of liquids have a high 
probability of responding to residual accelerations. Indeed, 
unexpected results have been attributed to residual accelera- 
tions in several cases [40, 86-88]. Particular examples of such 
experiments involve liquid (float) zones. Three types ot liq- 
uid zones are typically the subject of low-gravity experiments 
[89]: isothermal and non-isothermal zones suspended be- 
tween inert solids (i.e. liquid bridges), and non-isothermal 
zones suspended between a feed rod and a growing crystal. 
The dimensionless group associated with the shape of iso- 
thermal static liquid zones is the Bond number. Bo , defined 
earlier. It expresses the relative importance of gravitational 
and surface forces. A comprehensive description of stable 
zone shapes in zero gravity, and as a function of Bo , is given 
by Martinez etal. [41]. Fig. 10 displays their results for the sta- 



Fig. 10. Stability limits of an isothermal liquid bridge of volume V and 
length L suspended between discs of radius R and subject to a steady 
axial acceleration defined by a Bond number Bo, after Martinez etal 
[41], The angle 0 is defined in the inset 


Log B 



Fig. It. The maximum stable zone length of an isothermal liquid zone 
suspended between inert solids of radius R as a function of Bond num- 
ber: Comparison of theory * and experiment. After Coriell etal. [90] 


bility of static zones with a steady acceleration parallel to the 
zone axis. In order to relate the Bond number to a given ac- 
celeration level they give the value of Bo = 1.4 for a 11 vol- 
ume of water subject to an acceleration of 1 cm s~ 2 . For 
Bo = 0, the zone volume V — V* = n R 2 L and contact 
angle a = 90° the maximum stable length is 2n R (the Ray- 
leigh limit). Coriell etal. [90] investigated the effect of Bo on 
the maximum length of such zones and compared theoretical 
predictions with experimental results (see fig. 11). The actual 
length of the equilibrium zone is determined by the volume, 
the contact angle a and the static Bond number. 

Estimates of the effect on the free surface shape of oscilla- 
tory residual acceleration parallel to the axis of an isothermal 
liquid bridge have been made by Langbein [55]. Not supri- 
singly, his calculations suggest that the sensitivity of the sur- 
face is highest for disturbances with frequencies close to the 
natural frequency of the bridge, typically in the 0.001-10 Hz 
range. For the example given in [55], the maximum accelera- 
tion magnitudes that can be tolerated by a bridge less than 
90% of the maximum stable bridge length (given by the Ray- 
leigh limit) range from 5 - 10 _R to 5 * 1 0 ~ 3 g. The criterion 
used to determine the tolerable residual acceleration was that 
the amplitudes of all surface shape excitations with frequen- 
cies other than the resonance frequency be kept to 0.1% of 
the radius of the column. Examples of the results of this 
analysis are shown in fig. 12. The curves represent the maxi- 
mum tolerable axial acceleration for a given aspect ratio (fig. 
12 a), and for a surface response with a given number of 
nodes (fig. 12b). The results express the maximum tolerable 
frequency a Q in terms of the g-jitter frequency co , the reso- 
nant frequency, o) m , of the liquid zone associated with the 
m th node of the deformation of the zone surface, the viscosity 
of the liquid rand the spatial wavenumber 

a n < min ^ fco 2 - co}],) - a) r q ,„ j , (9) 

where or m = y q 2 m [(q m R) 2 - \}/qR, and 
q m - (m+ 1) n/H. 
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Fig. 12. Tolerable residual accferation fin units of g ** 
axial acceleration, following Langbein [55] 
a) Sensitivity for different aspect ratios A = L/'2 n R , 


R is the column radius, and 
H the column height. 

Note that the global minimum, and local minima of the 
tolerable acceleration curve occur at disturbance frequencies 
equal to the natural frequencies of the system. 

Zhang and Alexander [9\] have also analyzed this problem 
using a slender-body approximation. The problem of deter- 
mining the axisymmetric response of the shape of the free 
surface of a cylindrical liquid column bounded by two solid 
regions is modelled by a 1-D system of non-linear equations. 
It is found that the sensitivity of the zone to breakage and 



Fig. 13. Curves of tolerable acceleration fin units of g — 9,8 m s~ 2 ) vs. 
frequency for a maximum liquid bridge shape change of 10 °q at Bo * 
0.002, g n - t.42 • U)~ 5 g and C = 0.001. The solid curve , the dotted 
curve, and the dotted-dashed curve are (he results for aspect ratios L 2 
R = 2.6, 2.826. and 3.024, respectively . From Zhang and Alexander 
192] 
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shape changes depends on the static Bond number, aspect ra-_ 
tio and viscosity as well as the amplitude and frequency of 
the disturbance. The general trend is an increase in tolerable 
residual acceleration with increasing frequency. At the eigen- 
frequencies of the zone, however, there are strong deviations 
from this trend. At these frequencies this model also shows 
that the tolerable residual gravity level can be two orders of 
magnitude lower. These most sensitive frequencies have been 
found in the 10~ : — 10* 1 Hz range (see for example fig. 13). 
For the cases examined, maximum tolerable residual gravity 
levels as low as 10 -6 g have been calculated. For higher vis- 
cosities the tolerable acceleration level is increased for all fre- 
quencies. The equilibrium shape, as determined by the 
steady background acceleration, has a prounounced effect 
on the response at low frequencies. A change in slenderness 
of the bridge markedly changes the sensitivity to residual ac- 
celeration as the Rayleigh limit is approached. It is clear that 
for liquid bridges, small magnitude accelerations cannot be 
ignored. 

For non-isothermal zones, the dimensionless groups intro- 
duced earlier in sec. 3.1 govern the problem. For molten sili- 
con Sekerka and Coriell [89] have calculated the values of 
these groups. They are presented in table 4. The relative mag- 
nitudes of the terms suggest that on earth the zone shape will 
be dominated by the capillary effect (large Bo) for typical 
zone sizes, while under low-gravity conditions the zone 
shape will be influenced mainly by the temperature gradient 
along the interface. In addition (as discussed earlier) dynam- 
ic distortions may also be important. Considerations of the 
iloat zone crystal growth process introduce further complica- 
tions which ultimately require independent analysis. In the 
previous cases, the surfaces of the inert solids between which 
the zone is suspended are usually planar. This is generally 
not the case for the crystal growth process. The actual shape 
depends primarily on the thermal conditions. The influence 
of the meniscus angle at either of the crystal-melt-vapor tri- 
junctions can also be important and has been analyzed by 
Heywang [93] and Coriell etal. [94]. In addition to zone fail- 
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Table 4. Values of dimensionless groups for silicon calculated from the data of Chang and Wilcox [92]. The length scale is taken to be the zone radi- 
us R; If is estimated at 1.6- 10 K ~ 7 . For space, g k is taken as 10~ 4 g. At is assumed to be 10 C (after Sekerka and Conell [89]) 
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ure due to the classical capillary instability (or Rayleigh in- 
stability), a “dewetting” instability is found to be important. 
The instability arises when the angle between the meniscus 
and the melt-crystal interface is less than the wetting angle 
for the liquid melt on the solid crystal. For Si zones of greater 
than 1 cm radius, zone failure on earth occurs due to the 
Heywang instability [94, 95]. Zone failure by this mechanism 
is unlikely under low-gravity conditions. 

Feuerbacher et al. [56] developed an ad hoc model to de- 
scribe experiment sensitivity to residual accelerations. The 
model is based on the motion of a body in response to an ap- 
plied sinusoidal oscillation. The governing equation is that of 
a damped harmonic oscillator. The tolerable residual acceler- 
ation are (as expected) found to increase with increasing fre- 
quency of the disturbance. Furthermore, there is a minimum 
tolerable acceleration at disturbance frequencies correspond- 
ing to the natural frequency of the system. 

The effect of time dependent accelerations on experi- 
ments involving drops and bubbles has been examined for 
the cases of sinusoidal and Heaviside-function disturbances 
[57]. The criterion used to define the sensitivity is that the dis- 
placement of the drop not exceed 10% of the drop radius. 
For a given set of physical properties, the admissible acceler- 
ations can be expressed as a function of frequency (fig. 14). 



Fig. 14. Threshold intensity g * = gj R*/v : d T as a function of normal- 
ized frequency S = (0 R : 2v after Dewandre [57]. Here g r is the maxi- 
mum tolerable acceleration, d T the admissible displacement and r is the 
ratio of the drop density to that of the surrounding fluid 


The tolerable acceleration levels exhibit trends similar to 
those depicted in figs. 4-6. 

The influence of the gravity gradient, atmospheric drag 
and spacecraft attitude motion on the motion of a spherical 
solid particle immersed in viscous fluid has also been consid- 
ered [5, 6]. The extent and nature of the influence of the grav- 
ity gradient and the atmospheric drag on the trajectory of the 
particle are, for a given value of the Stokes drag coefficient, 
shown to be dependent on the characteristic orbital attitude, 
the magnitude of the atmospheric drag accelerations and the 
distance of the particle from the mass center of the space- 
craft. 

Amin [59] has investigated the influence of g-jitter on heat 
transfer from an isothermal sphere maintained at a constant 
temperature greater than the ambient fluid temperature. The 
body force was taken to have a single frequency sinusoidal 
time dependence and a small amplitude. A significant result 
of the analysis is that the Reynolds stresses associated with 
the fluctuating How result in a steady streaming motion. This 
has implications for both heat and mass transport in oscillat- 
ing flows. For the low Pr fluids examined, it was concluded 
that buoyancy induced convection caused by high frequency 
g-jitter (with amplitudes 10 -2 g) will result in heat transfer 
conditions significantly different from pure conduction. 
However, it was found that in high Prandtl number fluids 
with small kinematic viscosities, low frequency g-jitter will 
influence on heat transfer. 

Analyses of heat transfer in (laminar) oscillating flows in 
cylinders and between parallel plates [27-33] have shown 
that at high frequencies the effective thermal diffusivity is 
proportional to 


while for low' frequencies 

v 7/ ~ (o zSx : Wo. (11) 

Here co is the circular frequency, 

Jlx is the cross-stream average displacement of a fluid ele- 
ment over half the period of the oscillation, and 
Wo = L{( 0 /v) ] 2 is the Wormsley number which represents 
the ratio of the viscous response time to the period of the os- 
cillation. 

For a given frequency the heat transfer process was found 
to be most efficient for Pr = n/ Wo 2 i.e. when the character- 
istic time, r k = L 2 / a: associated with heat diffusion is equal 
to half the period of the oscillation. Fig. 15 illustrates this re- 
sult for fluids with different diffusivities and L = 1 cm. 
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W 0 = L(co/v) 1/z 

Fig. 15. Plot of effective diffusivity, vs. Wormsely number Wo, after 
Kurzweg [31] 

3.3 Linear Stability and Gravity Modulations 

The stability of a heated fluid layer of infinite lateral extent is 
affected by gravity modulation and has received a limited 
amount of attention for the case of sinusoidally modulated 
gravity [74-77]. The static case, where the gravity vector is 
perpendicular to the layer and parallel to the temperature 
gradient, is stable unless the critical Rayleigh number is ex- 
ceeded. The effect of gravity modulation is to introduce sta- 
bility for given combinations of modulation amplitude and 
frequency for situations which in the absence of modulations 
would exhibit instability. In contrast, for the case where the 
gravity vector and the temperature gradient are anti-parallel 
(which is stable for a steady gravity vector), instability can 
occur due to gravity modulation. 

The work of Biringen and Peltier [77] focused on the be- 
havior of a 3-D Ravieigh-Benard system subject to spatial 
and temporal gravity modulation. They used a pseudo-spec- 
tral method to examine the system numerically and studied 
the effects of sinusoidal and random gravity modulations 
with the gravity vector oriented perpendicular to the fluid 
layer. It was found that when the base acceleration was re- 
duced from 1 g to zero, the system parameters that had re- 
sulted in a synchronous response at 1 g led to a conductive 
state for the modulated zero g case. The spatially random 
modulations about a zero g base acceleration led to excita- 
tions of the 3-D temperature field. These excitations were 
less pronounced than for the 1 g case. For temporally ran- 
dom modulations about zero g, the cases considered that 
where excited at 1 g no longer exhibited excitation. 

The effect of oscillating accelerations and impulse forcing 
of an interface separating two immiscible fluids has recently 
been considered [78, 79], Both these analyses are restricted to 
time-dependent accelerations with zero mean which are 
oriented perpendicular to the interface between the fluids. 
The linear stability of a basic state characterized by a planar 
interface and an oscillating pressure. For the case of an arbi- 
trary oscillating disturbance it is shown that the linear evolu- 
tion of an infinitesimal perturbation of the interface is go- 
verned by a single integro-differential equation. Except in the 
limit of zero surface tension the effect of viscosity is shown 
to be small. As surface tension becomes dominant, resonant 
instabilities are excited at lower and lower wavenumbers. 
The interface is found to be more unstable to impulse type 
forcing than to sinusoidal forcing [79]. 


3.4 Direct Numerical Simulation 

The effects of low-gravity on the transport of heat and mo- 
mentum have been examined in a number of articles [37, 
60-71]. Robertson etal. [60, 61] found that for convection in 
circular cylinders with azimuthal variations in temperature at 
the boundary and the gravity vector applied perpendicular to 
the cylinder axis, the intensity of convection follows the pre- 
diction of Weinbaum’s first order theory [96] for low Ray- 
leigh numbers. Weinbaum’s First order theory predicts a sim- 
ple sinusoidal dependence of the maximum velocity as a 
function of orientation of the gravity vector. The effects of a 
variety of acceleration vectors (sinusoidal, cycloidal and lin- 
ear periodic) on motion in three fluids corresponding to mer- 
cury, helium and water (10“ 2 < Pr < 10) were studied by 
Spradley etal. [62]. They found that the steady mean part of 
the applied disturbance is more important than the oscillato- 
ry part (frequency = 1 Hz) in determining the flow field and 
heat transfer. Kamotani etal. [37] solved a linearized approxi- 
mation of the Boussinesq equations and investigated the ef- 
fect of an applied acceleration consisting of a time mean part 
and an oscillatory part on the temperature and flow fields in 
a rectangular enclosure. They found that the thermal convec- 
tion was predominantly oscillatory in nature. In addition 
they found that the acceleration perpendicular to the temper- 
ature gradient is the most important for the generation of 
fluid motion. 

The response of natural convection of a Boussinesq fluid 
in a cylinder to residual acceleration has been examined by 
Heiss etal. [64], and Schneider and Straub [ 97]. In the work of 
Heiss etal. the calculations are three-dimensional and in- 
volve gravity pulses and a rotating gravity vector. The ends of 
the cylinder are held at different temperatures. The cylinder 
walls are taken to be either perfectly insulating or perfectly 
conducting. The results are presented in terms of the relevant 
dimensionless parameters in the original papers. Fig. 16 illus- 



Fo* 0.02 0.04 


Fourier Number 

Fig . 16. Dimensionless maximum velocity caused by impulse distur- 
bances as a function of Fourier number Fo = t tc/d, after Heiss etal. 
[64]. Here t is the time [s] K is the thermal diffusivity [cnf s~ J ], d is the 
diameter of the cylinder [cm]. The parameter Ra * is the Rayleigh num- 
ber associated with the disturbance which is of dimensionless duration 
F f} * = 5.5- J0~ \ Ra* = fir g L? AT/v k, where L is the cylinder 
length [cm], AT is the temperature difference between the ends of the 
cylinder [K], (f Is the coefficient of thermal expansion [K~ 7 - g IS the 
magnitude of the acceleration [cm s~ 2 ] and v is the kinematic viscosity 
[cnr s ~ ; / 
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trates their results for impulse disturbances. The physical 
conditions were taken to correspond to: a cylinder of aspect 
ratio l, with diameter = 1 cm, v = 5 ■ 1 0 - 3 cm 2 s"\ k = 
7.14- 10" 3 cm 2 s’* 1 , P = 2.5 • 10” 4 , and a temperature differ- 
ence, AT, of 25 K. The magnitude of the acceleration pulse 
(duration = 0.7 s) corresponding to a Rayleigh number, 
Ra *, of 2,000 is - 3 • 10” 3 g results in a maximum velocity of 
approximately 5.8 cm s“*. The decay time for the flow field 
is about 5.6 s. Notice that for all the results, the decay time of 
the pulse is independent of the magnitude of Ra * For a ro- 
tating gravity vector of a given magnitude, the maximum ve- 
locity of the system is independent of the dimensionless fre- 
quency / = /* d 2 /K , up to / — 10 whereupon it decreases 
rapidly with increasing frequency. For the fluid with the 
physical properties described in the previous example, /* 
corresponds to 7.14* 10 -2 s“ 1 . 

Schneider and Straub [97] examined the effect of Prandtl 
number on the response to g-jitter for the same system as 
Heiss etal. They examined the response of the system to 
pulses of infinite length and sinusoidally varying pulses. For 
sinusoidal pulses they found for systems with Pr = 0.71 and 
7.0 (air and water) that the maximum velocity decreased with 
increasing frequency for dimensionless frequencies, F — f 
1 ?/ k , greater than 10. At Pr = 7 a maximum in the velocity 
was observed between F — 10 and 100 for the maximum am- 
plitude pulse examined (Ra = 2- 10 5 , or \g j « 3- 10 -4 g for 
a fluid with the properties of water, a 1 cm diameter cylinder 
and a temperature difference of 10 K). They determined that 
the system was least sensitive when the transient accelera- 
tions acted along the cylinder axis. 

The extent to which gravity causes buoyancy-driven fluid 
motion (and thus, solute redistribution) during directional 
solidiflcaiton has been examined using numerical models of 
buoyancy-driven convection in cylindrical and rectangular 
geometries [9, 65-69, 84, 98-101]. Except for [9, 65-69] these 
studies are restricted to axisymmetric situations in which a 
steady gravity vector is oriented parallel to the axis ol a cylin- 
drical ampoule. McFadden and Cone!! [66] have undertaken 
2-D calculations of the effects of time-dependent accelera- 
tions on lateral compositional variations during directional 
solidification. The gravitational acceleration was assumed to 



Fig. 17 . Maximum lateral solute segregation Ac/c^ {where c x is the 
far-; field solute concentration in the melt) in the crystal as a function of 
the period of the gravitational oscillation. After McFadden and Coriell 
[661 


have a uniform magnitude and rotation rate. The amount of 
compositional non-uniformity was found to increase with de- 
creasing rate of rotation (see fig. 17). 

Two- and three-dimensional models of directional solidifi- 
cation from dilute gallium-doped germanium melts have 
been used to determine the sensitivity of crystal composition- 
al uniformity to both time-dependent and steady residual ac- 
celerations [9, 67, 69]. The specific boundary conditions, ther- 
mo-physical properties of the melt, growth rates and am- 
poule size are all found to play a role in the determination 
of the experiment sensitivity. For a given set of operating 
conditions, it is found that at growth rates on the order of 
6.5 mm s _I the orientation of the experiment with respect to 
the steady component of the residual gravity is a crucial fac- 
tor in determining the suitability of the spacecraft as a means 
to suppress or eliminate unwanted effects caused by buoyant 
fluid motion. The worst case appears to be when the acceler- 
ation vector is parallel to the crystal interface. At growth 
rates on the order of microns per second, this orientation 
leads to compositional non-uniformities of 10-20% when the 
magnitude of the acceleration is of the order 10 -6 g. If, how- 
ever, the growth rate is lowered by an order of magnitude, 
the non-uniformity is reduced significantly (down to 1-5% in 
this case). Table 4 summarizes the results obtained for steady 
accelerations. 

The directional solidification (Bridgman-Stockbarger) 
process is also extremely sensitive to transient disturbances. 
The response of the system to a variety of impulses will be 
discussed. For example, a 3 - 10 -3 g impulse of one second 
duration acting parallel to the interface of a growing crystal 
produces a response in the solute field which lasts for nearly 
2,000 s. Consequently lateral and longitudinal compositional 
variations occur over a length of nearly 6 mm in the grown 
crystal. 

The response of the solute field and the lateral non-uni- 
formity to oscillatory accelerations varies from no response 
at all (at frequencies above 1 Hz with amplitudes below 
10- 3 g) to a significant response at 10' 3 Hz at amplitudes on 
the order of 10 -6 g. In addition, additive effects were ob- 
served for combinations of a steady and a low frequency re- 
sidual acceleration component. These additive effects gave 
rise to significant lateral and longitudinal non-uniformities in 
concentration. A number of different types of periodic dis- 
turbances were employed. Single frequency disturbances of 
the form g(t) = go_Fg„ cos(2 nf n t) were examined withgb = 0, 
\fl • 10" 6 and y 2 • 10~ 5 g , oriented parallel, perpendicular 
and at 45° to the crystal-melt interface. The frequency range 
examined was/, = 10~ 4 , 10~ 3 , 10' 2 , 10“ 1 and 10 Hz. For 
frequencies greater than 10 -2 Hz, there were no discemable 
effects on the solute fields. The velocity field did, however, 
respond to the oscillatory disturbances. For the case of 10 -3 
Hz (at 5- 10~ 6 g) the response of the solute field was signifi- 
cant. Lateral and longitudinal non-uniformity levels in excess 
of 15% were calculated. Fig. 18-20 show the lateral non-uni- 
formity as a function of time and highlights the additive ef- 
fect of oscillatory and steady components of the residual ac- 
celeration. 

The effect of a multiple frequency disturbance is illustrat- 
ed in fig. 20. The acceleration consists of steady and periodic 
contributions with the form: g(t) = go 4- g\ cos(2rc 10”%) 
+ gj cos (2)i 10--’ 0- Here |g 0 j - V^ ' 10 ” 6 8, ki I = 

3 y 2 • 10“ 6 g and ]g 2 I - 3 \fl ■ lO -5 g. The magnitude of 
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Fig. 18. Lateral non-uniformity in composition, f plotted as a function 
of time for an oscillator y residual acceleration with a maximum magni- 
tude of 3 \/2 • t0~ 6 g and a frequency of J0~ ? Hz, acting parallel to (he 
cnstal-melt interface [67]. The initial state was purely diffusive 

Fig. 19. Lateral non-uniformity in composition, g, plotted as a function 
of time for a residual acceleration consisting of a steady part with mag- 
nitude \f~2 ■ 10~ 6 g and an oscillatory part with a maximum magnitude 
of 3 \ r 2- 10~ 6 g and a frequency of Hz. acting parallel to the 
crystal-melt interface [67], The calculation was started from a steady 


flow associated with a \f2 • 1 0 6 g acceleration acting parallel to the 
interface 

Fig. 20. Lateral non-uniformity in composition, c. plotted as a function 
of time for a multi-component disturbance consisting of a steady low g 
background plus two periodic components: g(f) = go -b gj cos (2 n 
10~ s t) + g 2 cos (2 7t l(T 2 t), where \ g 0 \ — ; Si ] “ 

j/J. l(r 6 , and \ g 2 1 «= i/J- HT 5 [67 j. The calculation was started 
from a steady flow associated with a -y[2- I0~ 6 g acceleration acting 
perpendicular to the interface 


the compositional non-uniformity ± is seen to vary with the 
frequencies of the acceleration. More recently, it has been 
shown [67] that, particularly at high frequencies, it is impor- 
tant to consider the long time behavior of the system, be- 
cause for high frequency disturbances ( [« 0.1-1 Hz) the tran- 
sient behavior of the system is more sensitive than the long 
time behavior. For example, the response to a 1 Hz 10~ 2 g 
acceleration was found to exhibit a concentration non-uni- 
formity on the order of 10% during the first 80 s whereupon 
it decayed slowly to 1 % of over a period of 2,500 s. In addi- 
tion, calculations have been made using sample acceleration 
time histories obtained on Spacelab 3 (SL-3). The results var- 
ied according to the type of SL-3 data that were input. For 


the example shown in fig. 21. The response of the solute Field 
was negligible, although the velocity field fluctuated on the 
order of the crystal growth rate. The system exhibited some 
response to transient disturbances (for example the struc- 
tural response to thruster Firings) that were superimposed on 
background low frequency acceleration. Low frequency 
(< 10“ 1 Hz) acceleration components with amplitudes of 
10“ 3 g or more were also found to give rise to undesirable 
non-uniformity. 

Griffin and Mokatef[ 72, 73] modelled the nature of melt 
convection in a Bridgman system subject to steady and un- 
steady axial acceleration. They also found that at conditions 
corresponding to low-gravity, the response time of the melt 



Fig. 21. Lateral non-uniformity I £ | as a func- 
tion of time for accelerations taken from a 
sample time series (see inset) constructed from 
data obtained on Spacelab 3 [69 j. The accel- 
eration consists of a repeated "noise ' segment 
{length 10 s, 10- j Hz < f < 1 Hz) and a 
thruster firing event (length 10 s, 10 7 Hz < f 
< 13 Hz). The latter is introduced at 10 and 
80s 
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Table 5. Compositional non-uniformity $ [%] for Ge.Ga [9,80]; £ = 
(Cwtax-CwjJ/Cw, where c s denotes the crystal composition at the melt- 
crystal interface 


residual 

acceleration 

magnitude 

[*] 

orientation 

N 

ampoule width [cm] 

1 0.5 2.0 

growth rate [mm s' 1 ] 

6.5 0.65 6.5 6.5 

A) 10- J 

t - 

80 

to - 5 

— 

92.7 11.9 12.0 


y 

70.9 11.3 


i 

6.4 0.95 

1 

o 

i 

3.2 


/ 

39 


4- 

54.2 

10“ 6 

4- i 

11.3 2.0 


/ 

8.0 


\ 

0.7 0.0 

B) 10- 5 

— 

22.6 64.5 

10- 6 

- 

2.3 


e< is the unit vector parallel to g , Ais the normal vector to the crystal, 
I = 1 cm for all cases. Values in parentheses indicate 3-D results. 
All calculations were undertaken using the thermo-physical proper- 
ties of gallium-doped germanium (9) and A) and B) refer to the oper- 
ating conditions listed below. 

Operating conditions 

A) hot zone temperature ( T h ) 1,331 K 

distance between inlet and interface (L) 1.0 cm 
height of adiabatic zone (L a ) 2.5 mm 
ampoule width (diameter) 1.0 cm 

B) hot zone temperature {T tl ) 1,251 K 

velocity is proportional to the momentum diffusion time (i.e. 
it is controlled by the characteristic system length scale and 
melt viscosity). In accord with the general trend exhibited by 
most physical systems discussed in this article, they found 
that the velocity response to sinusoidal g-jitter decreased as 
the inverse square of the momentum diffusion time. 

A general conclusion that can be drawn from all attempts 
to characterize gravity-driven convective effects on direction- 
al solidification from two component melts is that the maxi- 
mum lateral solute non-uniformity (radial segregation for the 
axisymmetric cases) occurs near the transition from diffusion 
dominated to convection dominated growth conditions [52, 
84, 9]: that is, when convective velocities are of the same or- 
der of magnitude as the diffusive velocities. The conditions 
under which this “transition” takes place will depend on the 
specific nature of the forces driving convection. The orienta- 
tion of the steady component of the gravity vector is crucial 
in determining the magnitude of the gravity vector at which 
this transition occurs. Thus, for a given set of operating con- 
ditions the orientation of the gravity vector determines the 
suitability of a low-gravity environment for directional solidi- 
fication experiments. 

Motivated by the growth of friglycine sulfate (TGS) crys- 
tals on Spacelab 3, [11], Nadarajah etal. [102] have examined 
the acceleration sensitivity of crystal growth from solution. 
Thermal and solutal convection were included in their analy- 
ses and the crystal growth rate was chosen as the sensitivity 
parameter for the response to convective transport. Simula- 
tions were carried out for steady, impulsive and periodic ac- 
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celerations in order to determine tolerable acceleration lev- 
els. Long-time simulations of the experiment were conducted 
with steady background accelerations of 10~ 6 and 10“ 5 g. 
Impulsive and periodic disturbances of higher magnitudes 
were imposed at intermediate points. For steady accelera- 
tions the system was shown to have diffusion dominated 
transport at 10 -6 g but for 10" 5 g a transition to convection 
dominated transport occurs. Fig. 22 summarizes the results 
for oscillatory disturbances. For each case, the numbers 
correspond to the ratio g P Ac/ f v*, which is an estimate of 
the ratio of buoyancy to inertial forces. Here ft is the solute 
expansion coefficent, Ac the characteristic concentration dif- 
ference, /is the frequency of the acceleration and v* is the 
calculated characteristic fluid velocity. Clearly, for this exam- 
ple, a simple estimate alone will not suffice to determine the 
system sensitivity. The system is relatively stable to impulsive 
and periodic disturbances unless their magnitudes are very 
large. 



Fig. 22. Sensitivity of crystal growth from solution to various periodic 
disturbances. The solid line separates tolerable responses (10% or less 
growth rate fluctuations ) from intolerable ones. For each case, the num- 
bers correspond to the ratio gfi A c/f v* f which is an estimate of the the 
ratio of bouyancy and inertial forces. After Nadarajah etal. [102] 

4 Summary and Discussion 

The work described in this paper examines the sensitivity of 
a variety of space experiments to residual accelerations. In 
all the cases discussed, the sensitivity is related to the dynam- 
ic response of a fluid. In some cases the sensitivity can be de- 
fined by the magnitude of the response of the velocity field. 
This response may involve motion of the fluid associated 
with internal density gradients, or the motion of a free liquid 
surface. For fluids with internal density gradients, the type of 
acceleration to which the experiment is sensitive will depend 
on whether buoyancy-driven convection must be small in 
comparison to other types of fluid motion (such as thermoca- 
pillary flow), or fluid motion must be suppressed or elimi- 
nated (such as in diffusion studies, or directional solidifica- 
tion experiments). In the latter case, the experiments (for 
example diffusion experiments and directional solidification 
experiments) are sensitive to steady and low frequency accel- 
erations (< 10“ : Hz). For experiments such as the direction- 
al solidification of melts with two or more components, de- 
termination of the velocity response alone is insufficient to 
asses the sensitivity. The effect of the velocity on the compo- 
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sition and temperature field must be considered, particularly 
in the vicinity of the melt-crystal interface. As far as the re- 
sponse to transient disturbances is concerned the sensitivity 
is determined by both the magnitude and frequency the ac- 
celeration and the characteristic momentum and solute diffu- 
sion times. 

The directional solidification of doped germanium crys- 
tals has been carried out on Skylab [17, 18, 103, T04] and the 
Apollo-Soyuz Test Project (ASTP) [105, 106]. The latter in- 
volved the directional solidification of seeded Ge melts. 
Melts were doped with gallium and some were doped with 
1 % Si and 0.001 % Sb [105]. The results of these experiments 
revealed strong asymmetric non-uniformities in the space 
grown crystals. Lateral variations were also observed in sam- 
ples grown under terrestrial conditions but were much less 
pronounced [25]. It has been argued that the asymmetric re- 
distribution of the dopant can be ascribed to “barometric dif- 
fusion'’ of the solutes due to the acceleration gradient in the 
melt arising from the rotational motion of the spacecraft 
[105, 106]. The basis of the argument, however, appears to ig- 
nore the presence of gravity gradient and atmospheric drag 
effects and does not explicitly account for the spacecraft atti- 
tude motions. It has been demonstrated [9] that melt convec- 
tion can occur in response to a low magnitude steady residu- 
al acceleration, and that whenever the acceleration vector is 
not aligned with the ampoule axis strong asymmetries in 
composition can occur. Since there is a non-linear depen- 
dence of compositional uniformity on convection, there is al- 
ways the possibility that the observed non-uniformities in the 
terrestrially grown samples will be smaller than those ob- 
served in space grown samples. Thus, convection caused by 
residual acceleration cannot be discounted as the origin of 
the non-uniformities. 

Thermocapillary How experiments require that buoyancy 
effects are small in comparison to capillary effects. The esti- 
mated range of accelerations to which they are sensitive lies 
between 5* \0~ } Hz and 2 Hz for the most sensitive cases. 
Unless the magnitude of the effective buoyancy force (ex- 
pressed through the Grashof number) associated with a given 
acceleration is on the order of or greater than Rs (1 4- 
S, : ) I/2 , buoyancy effects are not likely to be significant. It is 
possible that more subtle effects of the acceleration environ- 
ment might occur (for example influencing the transition 
from steady to oscillatory flow) but this cannot be deter- 
mined by order of magnitude analysis. It is more likely that 
thermocapillary How experiments will be most sensitive to 
the motion of the free surface in response to g-jitter. 

The response of free surfaces to g-jitter has been ex- 
amined through analytical and numerical models. For cylin- 
drical liquid bridges this has been undertaken only for axial 
accelerations. The shape of the surface is most responsive to 
disturbances close to the natural frequency of the bridge and 
its harmonics. The sensitive frequencies lie in the 5* 10 _: * - 
10 Hz range, with maximum tolerable magnitudes from 
10“ 8 - 10“ : ' g depending on the sensitivity criterion and the 
particular model employed. 

Fig. 23 depicts the range of measured acceleration ampli- 
tudes (courtesy of Dr. H. Hamacher, DLR, and A/. J. B. Rog- 
ers, University of Alabama in Huntsville) as a function of fre- 
quency for several Shuffle missions. It should be remembered 
that in addition to the accelerations shown here there are 
gravity gradient and atmospheric drag accelerations with 



a Sled 

* Hop 

* Drop 

° Treadmill 
A Quiet 
a FPM ops 

* Stowage 
A PRCS 

0 V RCS 
O Drag 
v Rotation 
■ SL 3 < I Hz 


Fig. 23. Residual accelerations measured on orbit fin units of g — 
9,8 m s~ : j. f Data courtesy of Dr. H. Hamacher. DLR . and M. J. B. 
Rogers Univ. of Alabama in Huntsville). The Spaeelab 3 (SL-3) data 
are restricted to measured frequencies between JO'- and 1 Hz. The 
points corresponding to Sled, Hop , Drop, Treadmill , Quiet . FPM ops . 
Stowage, refer to activities and experiments on the D-J mission. FPM 
ops stands for Fluid Physics Module operations, P RCS and V RCS 
refer to the priman’ and vernier thrusters. The drag and rotation entries 
correspond to accelerations arising from slow variations in atmospheric 
drag during an orbit and attitude changes involving rotation 

magnitudes on the order of 10 -6 - 10 -7 g which can be 
steady (for gravity gradient stabilized type attitudes) or have 
frequencies that are twice the orbital frequencies (for a solar 
inertial attitude). In fig. 24 we depict these together with g - 
tolerance curves selected from the examples discussed ear- 
lier. There is some overlap between the sensitivity curves and 
measured accelerations. Even if the fact that most of the 
curves are obtained for order of magnitude estimates is taken 
into account, it is clear that a careful evaluation of the inter- 
action between experiments and the residual acceleration en- 
vironment is necessary. One the experiment sensitivity has 
been determined (which may sometimes require a more de- 



Log Frequency [Hz] 

Fig. 24. Acceleration sensitivity curves from selected examples discussed 
earlier superimposed on fig. 23. Cunes / and 3 correspond to the cases 
n = J. A — 0.9999 and n = l, A = 0.9 in fig. 12a . curve 2 corre- 
sponds to n — 2 and A = 0.9999 in fig. 1 2 b. Curve 4 corresponds to 
the most sensitive case of zone shape change in fig. 13. and curves 5 
and 6 respectively correspond to the semi-conductor and metal melt 
growth experiments from fig. 4. Curve 7 represents the thermo-diffusion 
experiment of fig. 5 
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tailed analysis than any discussed in this review) it should be 
possible for the experimenter to choose operating conditions 
which will minimize the effect of the predicted acceleration 
environment. Since flight opportunities are limited, it is par- 
ticularly important that the chances of success for any given 
experiment be greatly improved by detailed modelling prior 
to flight. This will allow investigators to better anticipate un- 
desirable effects of the prevailing residual acceleration con- 
ditions, and will ultimately lead to a better overall under- 
standing of the physical system. 
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Sensitivity of liquid bridges subject to axial residual acceleration 
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It has become evident that the float zone crystal growth method and isothermal liquid bridges 
may be very sensitive to the residual acceleration environment of a spacelab. Indeed, the shape 
equilibria show a high degree of sensitivity and, thus, even the small steady acceleration 
associated with the effective low gravity environment of a spacecraft cannot be ignored. Using 
a slender-body approximation, the problem of determining the axisymmetric response of the 
shape of the free surface of a cylindrical liquid column bounded by two solid regions is 
modeled by a 1-D system of nonlinear equations. It is found that the sensitivity of the zone 
shape depends on the static Bond number, B 0i aspect ratio, and viscosity, as well as the 
amplitude and frequency of the disturbance. The general trend is an increase in tolerable 
residual gravity with increasing frequency. At the eigenfrequencies of the zone, however, there 
are dramatic deviations from this trend. At these frequencies the tolerable residual gravity level 
can be two orders of magnitude lower at this frequency. For the cases considered the values of 
B 0 were taken to be 0.002 and 0.02 and the dimensionless viscosities, C = v(p/yR 0 ) 1/2 , used 
were in the range 0.001-0.01. Aspect ratios ranging from 96.3% to 82.7% of the Rayleigh 
limit were examined. For these cases, the frequencies associated with the lowest tolerable 
acceleration have been found in the 10“ 2 -10 _I Hz range. In terms of previously recorded and 
predicted residual accelerations, the sensitive frequency ranges for the cases examined are 
10~ 2 -10 _ 1 Hz and 1-10 Hz. Maximum tolerable residual gravity levels as low as 10~ 6 g have 
been calculated. The effect of viscosity is seen to increase the tolerable acceleration level for all 
frequencies. The equilibrium shape, as determined by the steady background acceleration, has 
a pronounced effect at low frequencies. A change in slenderness of the bridge markedly 
changes the sensitivity to residual acceleration as the Rayleigh limit is approached. 


1. INTRODUCTION 

The low gravity environment of a spacecraft affords the 
opportunity to eliminate or minimize the often undesireable 
effects of gravity on fluid behavior. Both the fluid physics 
and the materials preparation community have an interest in 

the sensitivity of liquid zones to residual accelerations. A 

truly zero-gravity environment is never achieved and the re- 

sidual acceleration that arises due to gravity gradient, atmo- 

spheric drag, crew motion, structural vibration, etc., 1 ' 3 is 

characterized by an ever changing magnitude and orienta- 

tion and affects the relative motion of any mass in the space- 
craft reference frame. Whether a particular space flight ex- 
periment involves the study of isothermal zone shapes, or 
attempts to utilize the low gravity environment to maximize 
the length of the liquid zone in a float-zone crystal growth 
experiment, the effects of residual acceleration must be con- 
sidered. 4 Experiments involving free surfaces of low viscos- 
ity liquids are sensitive to acceleration. Indeed, for some 
space experiments, unexpected results have been attributed 
to residual accelerations in several cases. 5 9 Martinez 5 and 
Haynes 6 comment on the effects of an attitude change dur- 
ing the study of liquid column stability and liquid spreading 
kinetics on Spacelab- 1. The residual acceleration environ- 
ment also had a noticeable effect on stationary fluid masses. 
Roll rates of 0. 13°/sec were apparently “...sufficient to inter- 
fere strongly with some of the fluid configurations stud- 
ied... .” 8 Later experiments on D-l that involved static and 


rotating liquid columns also experienced some interference 
that has been attributed to the residual acceleration environ- 
ment. 7 In particular, Martinez recorded the following inter- 
ference apparently caused by the residual acceleration envi- 
ronment. During the initial stages of the establishment of a 
long cylindrical column by liquid injection, a large ampli- 
tude deformation of the bridge was observed in response to 
shuttle maneuvers. Following this, a 10 cm liquid column 
exhibited a pronounced random oscillation. 

In an experiment seeking to excite a c-mode deforma- 
tion 7 of the bridge an amphora mode 7 resulted. The residual 
acceleration environment seems to have caused a shift in the 
stability limits. The maximum residual acceleration record- 
ed during these experiments was 10“ 4 g. 

In Table I we have expressed the Bond number 
B n = ( pgR l/y)g n /g> where y is the surface tension, R 0 is 
the radius of the ends of the zone, p is the liquid density, g n is 
the acceleration magnitude, and g — 980 cm sec 2 , for the 
physical properties of liquid columns used in two of the 
space experiments discussed earlier. These estimates indi- 
cate that the experiment involving the long liquid bridges 
should be more sensitive. In contrast to Martinez, none of 
Padday’s experiments were reported to have been affected by 
the residual accelerations. 

The numerical modeling of the response of a free liquid 
surface to time-dependent residual accelerations is compli- 
cated by the nonlinearity introduced by the free boundary 
motion. Few studies have been conducted on this subject. 
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TABLE I. Dynamic Bond numbers estimated for two D-l liquid zone ex- 
periments in which 5 cS silicone oil was used. 



Martinez 5 

Padday 6 

p(gcm~ 3 ) 

0.92 

0.92 

L(cm) 

10 

1.5 

pCdyn cm ') 

0.2 

0.19 

Bn 

4.5(10 ) 4 g n /g 

7.1(10 )~ i g n /g 


Furthermore, the models to date are based on the assump- 
tion of an inviscid bridge, 10 - 11 or have dealt with a linearized 
1-D model. 12 Thus, we have formulated a simple nonlinear 

I- D model of a viscous isothermal liquid zone and used it to 
examine the effect of g jitter oriented parallel to the zone 
axis. This allows an examination of interaction between var- 
ious components of a typically complex residual acceleration 
without the restriction that the response necessarily be lin- 
ear. The model is based on a formulation developed by 
Lee, 13 which has been successfully applied to several types of 
liquid bridge problems. 10,11,14-18 

We have examined the effects of single frequency ^-jitter 
disturbances. The zone response is described in terms of its 
deviation from its equilibrium shape. The equilibrium shape 
depends on the value of the steady background acceleration. 
The sensitivity is expressed in two ways: first zone breakage , 
defined in terms of the conditions under which the zone radi- 
us goes to zero, and second in terms of a shape change crite- 
rion , which indicates when the zone deviates more than 5% 
from its equilibrium shape. An optimal searching scheme is 
used to find the sensitivity limits of the bridge. 

II- FORMULATION OF THE MODEL EQUATIONS 

Consider a slender liquid bridge (see Fig. 1) held 
between two parallel rigid circular disks of radius R 0 sepa- 



FIG. i. Idealized liquid zone configuration. 


rated by a distance L. The disks are aligned coaxially. The 
liquid is an isothermal incompressible Newtonian fluid with 
constant physical properties. It is held in place by capillary 
forces. Consistent with our aim to isolate the mechanical 
response of the surface to an acceleration g(t) acting along 
the cylinder axis we make the following additional assump- 
tions: 

( i ) Internal motion of the liquid bridge is caused only by 
capillary pressure gradients caused by deformation of the 
surface. 

(ii) The effect of the atmosphere around the bridge is 
negligible. 

(iii) The interaction of the residual acceleration envi- 
ronment with the experiment apparatus is such that the top 
and bottom disks vibrate in phase and the frame of reference 
for the calculations is then conveniently taken to be attached 
to the disks. 

The governing equations are rendered nondimensional 
using R 0t (pR l/y) U2 , and R 0 (pR \/y) “ 1/2 to, respectively, 
scale distance, time, and velocity. Here p is the mass density 
and y is the surface tension. The dimensionless equations for 
mass and momentum transfer then take the following form. 

In the liquid (0 <r<R(z,t) f — A<z<A), 


1 d{ru) ^ dw 
r dr dz 


du du du 

b u b w 

dt dr dz 


dp * ^ d 2 u j 1 du d 2 u 

dr Wr 2 r dr dz 2 r V 


dw dw dw 

b U b w 

dt dr dz 


dp * r ( d 2 w 1 dw d 2 w \ 

dz \ dr 2 r dr dz 1 J ’ 


(1) 


( 2 ) 


(3) 


where p * is a reduced pressure given by p* = p -b B^{t)z , 
and C= v(p/yR) w 2 (with r the kinematic viscosity) is a 
measure of the relative strength of viscous and capillary 
forces. The function g(t) is given by 
g(t) = 1 4- (2?i/2? 0 )sin tot, where co is the dimensionless cir- 
cular frequency. Here, B 0 =z g^pR\/y is the usual static 
Bond number and B x = g^pR \/y is a dynamic Bond num- 
ber. We define the aspect ratio as A = £ /2R 0 , where L is the 
zone length. The magnitudes of the steady and time-depen- 
dent acceleration components are given by g 0 and g lf respec- 
tively. The linear frequency (Hz) of the disturbance is given 
by/= co(pR l/y) ’ 1/2 /2tt. 

The boundary conditions at the surface r = R(z,t) are 


P* 


2 C 

(1 +RI) 



1 

( 1+* 2 ) V2 


\+R] \ 

— g(t)z, 


(4) 


\ dr dz J \ dr dz > 


which, respectively, represent the balance of normal and tan- 
gential force components at the surface. In addition, there is 
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the kinematic condition 


dR , dR _ 
dt dZ 


(6) 


The conditions at z = ± A are 

/?(z,0 = 1, w(r,z,t) = u(r,z,t) = 0, (7) 

which fixes the contact line between the disk and the liquid 
and ensures no slip of the fluid in contact with the disk. At 
r = 0 the conditions 


w(0,z,0 = 0, -^-(0,z,0 = 0, (B) 

dr 

ensure axisymmetry. The formulation is complete upon spe- 
cifying suitable initial conditions that we take in the form 

R{z, 0) = R 0 (z), w{r,z,0) =u(r,z, 0) =0, (9) 

along with the constraint that the volume remain constant, 
i.e., 

J R 2 0 {z)dz = J R 2 (z,t)dz = 2\. (10) 

The above system of equations defines an unsteady free 
boundary problem since the location (and hence the shape) 
of the oscillating surface is a priori unknown and must be 
determined as part of the solution. 

In order to calculate the response of the shape of the 
liquid zone to time-dependent axial accelerations, we ap- 
proximate the system of equations presented in the last sec- 
tion and formulate a one-dimensional model. This consider- 
ably reduces the complexity of the problem and allows an 
examination of a relatively wide range of parameters without 
excessive computation times. Meseguer 10 used a one-dimen- 
sional model that is valid provided the slenderness A is large. 
(Recall that the slenderness of a static bridge is bounded 
from above by the Rayleigh limit, i.e., A max = 1 r.) If radial 
momentum effects are neglected 15 then Eqs. (2) and (3) 
become decoupled and the following system of equations 
results: 


1 d(ru) , dw_ _ n 
r dr dz 


dw , dw 

f- w 

dt dz 


dp* +c d 2 w 


(ID 


( 12 ) 


dz dz 2 

Since w is now depends only on z, Eq. ( 1 1 ) is easily integrat- 
ed and yields 

rdw 


u{r,z t t) = - 


2 dz 


(13) 


Following Lee 13 and Meseguer 10 we take 5 — R 2 (z,t) 
and Q = S(z f t)w(z i t) and substitute ( 13) into the kinematic 
boundary condition (6). This yields 


dt dz 


(14) 


which expresses conservation of mass at the surface. Similar- 
ly, Eq. ( 12) can be recast in the form 


3Q 

dt 


+ 


-(^) = 

dz\ s ) 


-S^- + SC—(^- 


')■ 


(15) 


dz dz 2 \S j 

which expresses the conservation of momentum in the axial 
direction for each cross section of the zone. At the surface, 


the component of force balance normal to the surface re- 
duces to 


p* = 


45+ ( dS/dz ) 2 

+ is^m 

dz dz l \S ) I 


(Kfr+KD 
(f)T 


3/2 


+ 4 45+i 


X 


25 + 




(16) 


where Eq. (14) and the expressions for 5 and Q have been 
used to replace u and w. Together with the following bound- 
ary and initial conditions, Eqs. (15) and (16) complete our 
approximate description of the physical system 


5( + A,t) = 1, Q{± A,/) =0, 

5( + A,0) =S 0 (z), G(±A,0)=0. 


111. SOLUTION METHOD 

The time-dependent response of the shape of a liquid 
zone to axial residual accelerations was calculated and 
curves describing the tolerable acceleration as a function of 
the frequency of the axial acceleration were obtained for a 
variety of conditions. The initial shape was determined from 
computed static shapes corresponding to the static Bond 
number B 0 using a method similar to that employed by Me- 
seguer and Sanz. 15 


A. Iteration procedure 

In order to facilitate the solution of the system of equa- 
tions (14)- (17) it is convenient to recast (15) and (16) in 
the form 


dt dz\ S ) 


- - 


p* = 


dz 

C 


+ CB{S) 


m), 

dz\S J 


45+ (dS/dz) 2 
+ 4[45 + 

X 


dSd 2 JQ\\ 
dz dzAs/J 

(f)T 


where 


->-(K(f)]M§)ri 


and 


B(S) 


HHf)'- 

*M§) 1 1 - 


dS cPS 
dz dz 2 


(18) 


(19) 


( 20 ) 


( 21 ) 


The solution scheme is based on a fully implicit time 
stepping procedure. However, since the position of the sur- 
face is unknown an iterative procedure is required to obtain 
solutions to the differential equations and the location of the 
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surface at each time. At the ( n + 1 )th time step the itera- 
tions proceed as follows. 

(i) Guess the surface shape for the initial iterate 

(flf = I). 

(ii) The pressure/?* at the {m + 1 )th iterate is calculat- 
ed in the following way. All terms in (19) are calculated 
from the previous (mth) iterate except for the leading vis- 
cous term. This is calculated from the solutions obtained at 
the previous [/2 th, (n — 1 )th] time steps using an Adams- 
Bashforth two-step approximation. 

(iii) Q is then obtained upon solution of the momentum 
equation (18) and boundary conditions (17) using a Briley- 
McDonald implicit method. 19 " 21 

(iv) The surface shape at the (m -f 1 ) th iterate is then 
obtained from a time-centered discretization of (14), 

(v) The surface shape is then compared with the initial 
guess (or the shape calculated at the mxh. iterate). The pro- 
cedure is repeated until the condition 

max(\Sr +] -S?\/S?)<E t (22) 

where is satisfied. Here a subscript “j” denotes the 
value of s at the /th point, while the superscript denotes the 
iterative step. The above procedure is repeated at each time 
step. 

The Briley-McDonald method used to solve (18) in- 
volves the evaluation of the nonlinear advection term at the 
(n - j- 1 ) th time step using 

^(g7S) y +1 = ^Q 2 /S) J + A; |<? ^(g7S) jj* 

+ 0[(Ar) 2 ] ( (23) 

where the time derivative has been approximated using a 
forward difference and the space derivative has been ap- 
proximated with a centered difference. A Crank-Nicolson 
scheme is used to increase the accuracy of the solution. The 
discretized form of (18) then becomes 

cr +l -g," / <?(g7S) V Af[ <?/ <?(£ 75) 

At \ dz ), 2 L<?zV dt )L 


1 

f(5^iy + 

■+ (s^yi 


2 

ll dz J, 

V dz J, J 


c 

„ d 2 (Q\ 

d/C 

)\] 

H 

A (S)~ (% 

1 + 5 S)-f(-f 

- 

2 

l dz 2 \S ) 

1 dz\S 

7 J 

c 

r „ a 2 /o\ 

i d(0 

>\1 

T - — 

A (S)^~ 

+ 5(S)-H-s 

- 

2 

dz 2 \S) 

dz\ S 

A 


where a centered difference is used for spatial discretization. 
Equation (24) is unconditionally stable and has second-or- 
der accuracy with a truncation error of O [ ( A*) 2 , (At) 2 ]. 

In addition to the boundary conditions (17) we also 
need to evaluate the pressure at the disks. This is achieved as 
follows. First we use the fact that ds/dt — 0 at z = ± A 
which implies dQ/dz = 0 [from ( 14) ] and that dQ/dt — 0 
atz = ± A. We can then obtain the pressure from (18) eval- 
uated at z = ± A. This requires the evaluation of dp/dz and 
d 2 Q/dz 2 at z = f A. The pressure derivative is approxi- 
mated using the following standard three-point forward dif- 
ference expression: 



Pi - 4 /h + fop 
± 2Az 


-f 0(Az 3 ), 


(25) 


where, p 0 =p(± A), p x = p{ ± A =p Az), 

Pi— Pi ± A + 2Az), and Az is the spatial step size. Similar- 
ly, at z = ± A, d 2 Q/dz z is evaluated using 


f d 2 Q \ 

\ dz 2 ) z - 


%Q\ — Qi — 7Qo 


-bO(Az 4 ), (26) 


± a 2Az 2 

where Q { = Q{±A + Az), Q 2 = Q( ± A + 2Az), and the 
fact that dQ/dz — 0 has been employed. Upon substitution 
of (25) and (26) into (18) and rearranging terms we obtain 
p 0 with second-order accuracy. 

The system of discretized equations and boundary con- 
ditions is solved using a Thomas algorithm (method of fac- 
torization). 19 


B. Initial iterate 

For this method of solution the choice of surface shape 
for the initial iterate at each time step affects the speed of 
convergence. For a poor choice the iteration may even fail to 
converge. To ensure that the initial iterate is reasonably close 
to the solution, the following formula is employed: 

<27> 

For the computations presented here the number of iter- 
ations required was six or less whenever (22) was satisfied 
for the initial iterate. 



FIG. 2. Computational flow chart. 
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TABLE II. Effect of changing Af on the tolerable acceleration g r 


O) 

A 

Bo 

gt 

0.001 25 

314 

2.6 

0.002 

0.255 

0.0005 

31.4 

2.6 

0.002 

0.256 


C. Optimal searching scheme 

In order to define the sensitivity limits for the liquid 
zone shape an optimal searching scheme was used to delin- 
eate the boundary between the regions in parameter space 
for which solutions either do or do not satisfy our sensitivity 
criterion. The scheme is described in the flow chart shown in 
Fig. 2. Eight to ten calculations were typically needed to 
obtain each point on the sensitivity curves. 

D. Accuracy of the numerical scheme 

An investigation into the accuracy of the results re- 
vealed that the number of time steps per period of the distur- 
bance should be at least 160. Increasing the number of steps 
to 400 per period did not, however, result in any significant 
change in our results even at high frequencies (see Table II ). 

The spatial accuracy also proved to be an important 
consideration. Results obtained using the MacDonald-Bri- 
ley scheme required 97 spatial points. 

In order to test the adequacy of the 1-D model we com- 
pared results from 1-D calculations to results obtained by 
Sanz 16 for breaking times of an inviscid liquid bridge subject 
to an initial deformation of the surface area S given by 

S(z, 0) = 1+6(6- 2)sin(77-z/A), Q(z, 0) - 0. 

The results of the comparison are given in Table III. Our 
calculations are in reasonable agreement with those of Sanz. 
Table IV lists a comparison between our results for the first 
eigenfrequency, those predicted by the linearized analysis of 
Meseguer, 1 1 and experimentally obtained values, 1 1 Our cal- 
culated values lie between the linearized results and the ex- 
perimentally obtained values. Finally, we note that the anal- 
ysis of Langbein, 12 employs an expression for a damped 
linear oscillator together with natural frequencies of infinite 
liquid columns obtained by Bauer 22 and will overestimate 
the natural frequencies of finite bridges. This has been well 
documented by Sanz 16 in a thorough comparison of various 
methods used to simulate liquid bridge phenomenon. 


TABLE III. Comparison of breaking times ( T b ) and partial breaking vol- 
ume ( V b ) with those of Sanz 1 6 (7+, F*). Here A is the slenderness, e is the 
initial perturbation amplitude, and Az is the spatial step size used in our 
calculations. 


A 

£ 

A z 

T b 

Tt 

V> 

vi 

3.07 

0.2 

0.128 

21.02 

21.98 

0.8525 

0.8554 

3.14 

0.2 

0.130 

13.71 

13.62 

0.8474 

0.8574 

3.20 

0.2 

0.134 

11.82 

11.52 

0.8427 

0.8442 

3.10 

0.3 

0.129 

7.97 

8.4 

0.8525 

0.8539 

3.10 

0.4 

0.129 

4.41 

4.84 

0.859 45 

0.8566 


TABLE IV. Comparison of the dimensionless first eigenfrequency co u com- 
puted using our method, the linearized method of Meseguer u {*)*), and 
experimental results" (&„,). 


A 

<y 0 

«8 

O) 0t . 

2.571 

0.312 

0.348 

0.30 

2.714 

0.248 

0.259 

0.23 

2.854 

0.174 

0.222 

0.19 


IV. RESULTS AND DISCUSSION 

The calculations were undertaken for a disk radius of 
0.0175 m and liquids with the following physical properties; 
p = 920 kg m“ \ y = 0.02 N m~ *, and kinematic viscosities 
in the range 6X 10~ 7 -6x 10 5 irr sec *. These viscosities 
cover a range of typical experimental materials (water, sili- 
cone oils, etc., used in experiments 10 ). The sensitivity crite- 
ria chosen to characterize the response of the bridge repre- 
sent two extremes. The first is the deviation of the the bridge 
shape by more than 5% from its equilibrium radius, the sec- 
ond is breakage of the bridge. The effects of viscosity, back- 
ground steady acceleration (static Bond number, B 0 ), and 
slenderness ( A ) were examined and the results are presented 
below. Each point on these curves was obtained after run- 
ning the calculation for times corresponding to more than 
ten periods of the driving force. 

Figures 3 and 4 show the effect of viscosity on the toler- 
able acceleration, g n for a frequency of 5 Hz at fixed A and 
B i)y for both sensitivity criteria. The tolerance increases as 
the viscosity is increased. Notice, however, that for values of 
the dimensionless viscosity parameter C less than 10 1 the 



FIG. 3. Dimensionless viscosity parameter C versus tolerable acceleration 
for breakage of the bridge at /= 5 Hz, A - 2.6, B t) ~ 0.002, and 
*,= 1.42X10'** 
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FIG. 4. Effect of viscosity on tolerable acceleration for the shape change 
criterion at/= 5 Hz, A = 2.6, B () = 0.002, and g () = 1.42 X 10 ' * g. 


FIG. 6. Curves of tolerable acceleration versus frequency for shape change 
criterion at B 0 = 0.002, g () = I.42X lO'* 5 g, and A = 2.6. The solid curve 
and the dashed curve are the results for C = 0.001 and 0.0 1, respectively. 


increase in tolerance is slight even though the viscosity has 
been increased by two orders of magnitude. Only when C 
approaches 10“ 1 does the increase in tolerance become sig- 
nificant. Figures 5 and 6 further illustrate this point. The 5 
Hz residual accelerations greater than 0.03 g would not nor- 
mally be expected to occur in spacecraft laboratories so we 
have not examined cases with C> 0. 1 . 


The value of the static Bond number has a more pro- 
nounced effect on sensitivity as seen in Figs. 7 and 8. This is 
limited, however, to the lower frequency range. With in- 
creasing frequency the differences are less pronounced. This 
is particularly so for the shape change criterion. It is interest- 



FIG. 5. Curves of tolerable acceleration versus frequency for breakage of 
the bridge at B {) = 0.002, g„ = 1.42 X 10' 5 g, and A = 2.6. The solid curve 
and the dashed curve are the results for C = 0.001 and 0.01, respectively. 



FIG. 7. Curves of tolerable acceleration versus frequency for shape change 
criterion at A = 2.6 and C = 0.001. The solid curve and the dashed curve 
are the results for B n = 0.002 (g 0 = 1.42X 10“ s g) and B {) = 0.02 
(g n — 1.42X 10“ 4 g), respectively. 


1971 


Phys. Fluids A, Vol. 2, No. 1 1, November 1 990 


Y. Zhang and J. 1. D. Alexander 


1971 








FIG. 8. Curves of tolerable acceleration versus frequency for breakage of 
the bridge at A = 2.6 and C = 0.00 1 . The solid curve and the dashed curve 
are the results for B ( , = 0.002 ( ftj = 1.42X I0“ 5 g) and B 0 = 0.02 
(g () = 1.42X 10~ 4 g), respectively. 


ing to note the additional deviations from the general trend 
for the case of zone breakage between 10 ~ 1 and 1 Hz. 

A small decrease in the slenderness A changes the mag- 
nitude of tolerable acceleration markedly when A is close to 



FIG. 9. Curves of tolerable acceleration versus frequency for breakage of 
the bridge at B 0 = 0.002, g 0 = .42 X 1 0 “ 5 g, and C = 0.001. The solid curve, 
the dotted curve, and the dot-dashed curve are the results for A = 2.6, 
2.826, and 3.024, respectively. 



FIG. 10. Curves of tolerable acceleration versus frequency for the shape 
change criterion at B 0 — 0.002, g 0 — 1 .42 X 1 0 “ 5 g, and C = 0.00 1 . The sol- 
id curve, the dotted curve, and the dot-dashed curve are the results for 
A = 2.6, 2.826, and 3.024, respectively. 


A max . This is illustrated in Figs. 9 and 10. The change in 
tolerable acceleration magnitude becomes less significant for 
values of A less than 80% of the critical value. 

Of most practical interest is the sensitivity of a given 
zone (i.e., A, B m and C fixed) as a function of the frequency 
of the disturbance. This is shown in Figs. 5-10. For both 
criteria the general trend is an increase in tolerable residual 
gravity with increasing frequency. At the eigenfrequencies 
of the zone there are, however, dramatic deviations from this 
trend. Here the liquid zone is extremely sensitive in compari- 
son to both higher and lower frequencies. The tolerable re- 
sidual gravity level can be two orders of magnitude lower at 
this frequency. For the values of A and B 0 (g 0 ) considered 
these frequencies are found in the neighborhood of 1 0 1 Hz. 
Associated maximum tolerable residual gravity levels as low 
as 10 -3 to 10~ 6 g have been calculated. In addition, at fre- 
quencies higher than the most sensitive one there are less 
dramatic deviations from the general trend. The frequencies 
at which these deviations occur are the eigenfrequencies for 


TABLE V. Most sensitive frequencies and associated tolerable acceleration 
(Hz) calculated using the 1-D (/ 0 , g to , ) model and Langbein’s if*, gf ol ) 
linear oscillator model. 12 For the 1-D model B 0 = 0.002 in ail cases. Lang- 
bein’s model is independent of B 0 . 


A /„( Hz) f* (Hz) &„,(*) (*) v(m 2 sec ') 

3.024 0.0334 0.0663 6.32X10- 7 1.11 X 10 9 5X10" 7 

2.826 0.0070 0.1224 1.16x10 ” 2.18X10" 9 5x10 7 

2.6 0.1049 0.1857 3.26xl0 _<> 3.6XI0" 9 5x10 

2.6 0.1049 0.1857 3.53x10 '’ 3.6x10"“ 5xl0" 6 


1972 


Phys. Fluids A, Vol. 2, No. 1 1 . November 1 990 


Y. Zhang and J. 1. D. Alexander 


1972 





TABLE VI. Tolerable acceleration (g , ) and associated eigenfrequencies (a?*) for two of the cases (shape change criterion). 


A = 3.024, B 0 — 0.002, C=0 
i o * 

g t (msec -2 ) 

0.21 

6 . 0 X 10- 6 

4.0 

6.09X10- 4 

11.0 

1 . 88 X 10" 3 

20.5 

5.10X10" 3 

33.5 

5.06X10- 3 

A = 2.826, B 0 = 0.002, C = 0 
g f (msec -2 ) 

0.44 

1.09X10 -5 

4.7 

7.79X10-" 

12.56 

1.96X10' 3 

24 

4.76X10- 3 

38.5 

5.78X10- 3 


the zone. As the frequency of the forcing function ap- 
proaches an eigenfrequency the zone becomes more sensi- 
tive. The lowest tolerable acceleration occurs in conjunction 
with the smallest eigenfrequency. We note that owing to our 
assumption that the top and bottom disks vibrate in phase, 
the eigenfrequencies locally exhibiting the most sensitive re- 
sponse are those associated with eigenmodes having an odd 
number of nodes (i.e. } an even number of half waves). This is 
in agreement with the analysis of Mesegeuer. 11 It is interest- 
ing to compare our results with linear oscillator model of 
Langbein. 12 Table V gives a comparison between the linear 
oscillator model and our 1-D model. While the results are 
qualitatively similar there are three distinct differences. 
First, for the case considered the linear oscillator model 
overestimates the sensitivities by almost two orders of mag- 
nitude. Second, as mentioned previously, the eigenfrequen- 
cies predicted by our model are smaller. The third difference 
is that Langbein’s model implicitly assumes that the axial 
vibration will excite deformations with a both odd and even 
number of nodes (this is precluded in our model). If the top 
and bottom disks do vibrate out of phase, then our results 
will underestimate the response of the liquid bridge at eigen- 



FIG. 11. Shape change as a function of time for a liquid zone with A — 2.6, 
C = 0.001, and g(t ) = 1 -j- (B { /B 0 ) sin (cot). The dimensional frequency, 
f is 0.06 Hz, and B t /B 0 = 12. For the physical properties considered, the 
maximum oscillating acceleration magnitude is 1.8x 10~ 4 g. 



FIG. 12. Shape change as a function of time for a liquid zone with A = 2.6, 
C— 0.001, andg(/) — l + (2?,/2? Q ) sin (cot). The dimensional frequency, 
/ is 3 Hz, and B t /B 0 = 916. For the physical properties considered, the 
maximum oscillating acceleration magnitude is 1 .3 X 10“ 2 g. 


frequencies corresponding to eigenmodes with an even num- 
ber of nodes. 

Table VI lists the eigenfrequencies (which can be deter- 
mined independently 11 ) and the associated maximum toler- 
able acceleration, g n for selected inviscid cases. For the in- 
viscid zones, we found that whenever the higher 
eigenfrequency is close to an integer multiple of a lower ei- 
genfrequency the sensitivity of the zone to disturbances with 
the higher frequency is increased. This noticeable increase in 
sensitivity is caused by nonlinear interaction and excitation 
of the lower mode. 

The spatial deformation of the liquid zone also varies as 
a function of the forcing frequency. Figures 1 1 and 12 illus- 
trate the evolution of the zone shape as a function of time for 
selected cases. Note that, owing to the small value of 2? 0 , for 
these cases the initial surface shape is visually indistinguish- 


° Sled 

• Hop 

• Drop 

° Treadmill 
A Quiet 
■ FPM ops 
+ Stowage 
A PRCS 
° V RCS 
O Drag 

V Rotation 

• SL 3 < 1 Hz 


-4 * 3 - 2-1012 

LOG FREQUENCY [Hz] 



FIG. 13. Residual accelerations measured on orbit (in units of g = 9.8 
m sec -2 ) and selected sensitivity curves, Curves 1 and 2 correspond to 
A = 2.826 and 3.024 from Fig. 10. Curve 3 corresponds to A = 3.024 from 
Fig. 9. (Acceleration data courtesy of Dr. H. Hamacher, DLR, and M. J. B. 
Rogers, Univ. of Alabama in Huntsville.) The Spacelab 3 (SL-3) data are 
restricted to measured frequencies between 10“ 1 and 1 Hz. The points cor- 
responding to Sled, Hop, Drop, Treadmill, Quiet, FPM ops (Fluid Physics 
Module operations), Stowage, refer to activities and experiments on the D 
1 mission. 2 ' P RCS and V RCS refer to the primary and vernier thrusters. 
The drag and rotation entries correspond to accelerations arising from slow 
variations in atmospheric drag during an orbit and attitude changes involv- 
ing rotation. 
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able from the surface of a right circular cylinder. While high- 
er forcing frequencies are associated with more nodes, the 
shape of the zone is complicated by the excitations of the 
fundamental mode and by nonlinear interactions. 

In summary, our results indicate that the zone is most 
sensitive to accelerations with frequencies close to or equal 
to the lowest natural frequency of the zone. It is useful to 
assess the sensitivity of the zone in terms of predicted space 
station and/or spacelab environments and with previous ac- 
celeration measurements. For the liquid zone cases exam- 
ined we have seen that frequencies in the 10“ 2 -10 -1 Hz 
range appear to be the most sensitive. Figure 13 shows sensi- 
tivity curves selected from Figs. 9 and 10 together with accel- 
eration as a function of frequency taken data recorded on 
orbit. The amplitudes of low-frequency ( < 1CT 2 Hz) acce- 
lerations predicted for the space station should not exceed 
levels of 10-V 3 Higher frequencies can be associated with 
acceleration magnitudes of up to 10 -2 g. 24 In terms of these 
predicted levels, or those measured on past missions, 3 25 the 
practical sensitivity range is restricted to disturbances with 
frequencies in the range 1 0 — 2 — 10 1 Hz and 1-10 Hz. 
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SURFACE TENSION AND BUOYANCY-DRIVEN FLOW IN 
A NON-ISOTHERMAL LIQUID BRIDGE 

YIQIANG ZHANG AND J. IWAN D. ALEXANDER 
Center for Microgravity and Materials Research , University of Alabama in Huntsville , Huntsville. AL 35899, U.S.A. 


SUMMARY 

The Navier-Stokes-Boussinesq equations governing the transport of momentum, mass and heat in a non- 
isothermal liquid bridge with a temperature-dependent surface tension are solved using a vorticity-strcam- 
funciion formulation together with a non-orthogonal co-ordinate transformation. The equations are 
discretized using a pseudo- unsteady semi-implicit finite difference scheme and are solved by the ADI 
method. A Picard-type iteration is adopted which consists of inner and outer iterative processes. The outer 
iteration is used to update the shape of the free surface. Two schemes have been used for the outer iteration; 
both use the force balance normal to the free surface as the distinguished boundary condition. The first 
scheme involves successive approximation by the direct solution of the distinguished boundary condition: 

The second scheme uses the artificial force imbalance between the fluid pressure, viscous and capillary forces 
at the free surface which arises when the boundary condition for force balance normal to the surface is not 
satisfied. This artificial imbalance is then used to change the surface shape until the distinguished boundary 
condition is satisfied. These schemes have been used to examine a variety of model liquid bridge situations 
including purely thermocapillary-driven flow situations and mixed thermocapillary- and bouyancy-driven 
flow, 

keywords 'ThermocofiUAry pow, buoyancy, free. Surface., f mht 

Picard iberaho*, 

1. INTRODUCTION 

The computation of solutions to the steady free boundary problem of mixed thermocapillary- and 
buoyancy-driven convection in a non-isothermal cylindrical liquid bridge is complicated by the 
strongly non-linear boundary conditions at the unknown free surface. Most studies to date have 
avoided the computation of the free surface shape and assumed that the liquid surface is a circular 
cylinder 1-6 or have imposed a non-circular cylindrical surface shape. 7 These models have 
involved either half-zone configurations (where the liquid bridge is held between rigid disks of 
different temperature) or full zones (where the liquid bridge is held between two solids of equal 
temperature). The full zone models are motivated by the floating zone crystal growth process. 8 
Recently, Duranceau and Brown 9 have approached the full zone crystal growth problem using 
the finite element method and have computed the shape of the liquid surface as well as the 
melt-crystal and melt-feed rod surfaces together with interacting thermocapiliary and buoyant 
convective flow. Lan and Kou 10 have also approached the full zone problem using a finite volume 
method but restricted their calculations to the zero-gravity case for which the surface deformation 
from the circular cylindrical shape is minimal and buoyancy-driven flow is absent. Hyer et al. 11 
have used a finite dement method, which is well suited to irregular geometries, to compute 
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interacting thermocapillary and buoyant convective flow in full and half-zone configurations bu 
did not consider solidification. 

Recently, finite difference methods have been used to solve problems with free and movinj 
boundary geometries using various mapping techniques. 12 These methods are also applicable tc 
the free boundary problem associated with liquid bridges and floating zones. Kang and Leal ' 
used the finite difference method with orthogonal boundary-fitted co-ordinates to study th( 
deformation of a bubble. The boundary-fitted mapping scheme requires the solution of a couplet 
set of Laplace equations to determine the new grid to be generated at each outer iteration. 

In the present investigation the governing equations are recast in terms of a vorticity-stream 
function formulation together with a non-orthogonal co-ordinate transformation. The latte 
allows an irregular free boundary to coincide with a co-ordinate line (or surface) without the neei 
to solve a coupled set of Laplace equations. The resulting equations are discretized usm; 
a pseudo-unsteady semi-implicit difference scheme and solved by the ADI method. The combina 
tion of the above methods provides a reasonably accurate and economical solution procedure 
Four boundary conditions are specified at the free surface: the kinematic boundary condition th. 
balance of energy across the surface and the balance of force normal and tangent to the surface 
The energy balance, tangential force balance and kinematic conditions at the free surface ar 
solved together with the Navier-Stokes and continuity equations, while the normal force balanc 
condition is distinguished 14 to determine the free surface shape. In addition, an ‘outer iteraliv 
procedure is needed to locate the free surface. In this paper two outer iterative schemes ar 
reported. The first scheme involves successive approximation by the direct solution of the lore 
balance normal to the free surface. The second scheme, after Ryskin and Leal, uses the fore 
imbalance between the fluid pressure, viscous and capillary forces at the free surface which anse 
when the boundary condition for force balance normal to the surface is not satis e . 1 

artificial imbalance is used to drive the surface shape towards its equilibrium position (re. unt 
the force balance condition is satisfied). These schemes are used to examine a variety of mod* 
liquid bridge situations including purely thermocapillary-driven flow situations and mixe- 
buoyancy-thermocapillary-driven flow. 


r \ 

j If k OKi Oi 
pinod 




rfef ’I 

2. FORMULATION OF THE PROBLEM 


2.1. Governing equations 

Consider a cylindrical liquid bridge (see Figure 1) held between two parallel coaxial circuit: 
rigid disks of radius R 0 separated by a distance L. The liquid is a non-isothermal incompressibi 
Newtonian fluid. The bridge is held between the disks by surface tension. The free surface of th 
bridge is a gas-liquid surface and the steady surface shape is described by r = R(z). Each disk 
maintained at a constant temperature T 0 . Surface heating is provided through an ambiei 
temperature r*(r). Radiative and convective heat transfer at the free surface are accounted for t 
a heat transfer coefficient h. In addition, we assume that the gravitational acceleration is paralk 
Wihl to the cylinder axis andflhe velocity a nd thnt ri w temperature field and the deformation of the fre 
A surface are axisymmetne. Furthermore, we let the surface tension at the free surface vary hnearl 
with temperature and assume that the Boussinesq approximation holds. 

The governing equations are made dimensionless by scaling length, time and velocity with R, 
R 0/ L'* and U* respectively. Here V* is a characteristic velocity given by 


ft 
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Figure 1. Liquid bndge model 


10 “ mp "“" re ““ ' “ tte 

dynamic viscosity. , , when th end disks are held at the maximum and 

We shall refer to a half-zone model whcn thc temperature maximum occurs 

minimum temperatures respectively and a full lempcratU rc r.« to be parabolic 

s — *-«*. - - 

The non-dimensional pressure is 

p* + Po9Z 0 

p ‘7^ n ' 

where P . , s „ — rr! -x: 

axial co-ordinate and p 0 »s the density P lhese sca les the dimensionless steady 


cT 


1 c(ru) + Sw =Q ^ 
r dr dz 

(1) 

CU fP.ip + i^ + S4, 
T z = ~7r ReW rdr dz 2 r* J 

(2) 

dz + Re \ dr 1 r dr dz 1 J Re 

(3) 

dT l (PT 

dz .VfaVtlr 2 r dr dz 2 J & 

(4) 
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where the Reynolds number Re ; Marangoni number Xfa and Grashof number Gr are respectively 




RqU* 


Xfa = 


!y|A77? 0 


Gr = 


gfiATRi 


v v 

Here v is the kinematic viscosity, k is the thermal diffusivity, fi is the volumetric thermal 
expansion coefficient and g is the gravitational acceleration. 


2.2. Boundary conditions 

At the disks the boundary conditions are 

u = vv=F=0 at : = ± A/2, 

and the symmetry conditions at the centreline r = 0 are 

_cHv_cT 
U dr cr 

The boundary conditions at the free surface r— R(z) take the form 


(5) 


( 6 ) 


2 Re 1 pu /i 

cR\ 2 dw dR j 


p-Gz + A- j [_ar + V 

dz J dz dz 


Re-'{Co l -T) | 

f 1+ (dR/dz ) 2 

c 2 R\ 

+ [1 +(<5f?/dz) 2 ] 3 ' 2 


dz 2 } ’ 


(7) 


1- 


' dR x2 
dz 


du dw\ dRfdu f, /g*VT /2 (¥L + \ (g) 

dz + ^) + 2 dz\dr dz) [ V^/ J U* dz dr / ’ 


dR A 
u + w — =0, 
cz 


— — ^) + Bi(T- r x )=o, 
dR/dz) 2 )' 1 ' 2 «/ 


(9) 

( 10 ) 


where 


(1 +(dR/d 

;AT 


C 0 = - 


B/= 




G = 


Vo 


l/* 1 


are the capillary number, Biot number and dimensionless gravitational acceleration respectively 
and 7o is the mean surface tension. The force balances at the free surface in the normal and 
tangential directions are given by equations (7) and (8) respectively. Equation (9) is t e lnematic 
boundary condition at the liquid- gas surface . The thermal boundary condition at the surface is 
given by equation (10) with th^ a qwTvai ^ heat transfer coefficient Bi. The constant A on the 
left-hand side of (7) represe;u»^dimensionless reference pressure difference across the surface, n 
liquid bridge modeUysfems with fixed rigid disks such as the one discussed here, A is determined 
by the cgD£tfltff~volume constraint 

y 0 = j" 7 rR 2 (z)dz = constant. 


mensionless 


j - A/2 



Finally, the condition that the contact 
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lines between the liquid end disks are fixed is 


R = 1 at :=±A/2. 
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( 12 ) 


3. numerical method 

J _1 r)\ii 


l C\jj 


U~- -7-1 

r dz 


(D = 


cu 

dz 


w=- 

dw 

^'dr' 


— 1 d\ft 

Hr' 


From (13), (14), (2) and (3) one obtains for oj 


COJ c 
u— + vv , 
cr dz 


’CO cou 1 /*’« lto ±« * 

: = TeW rcr + dz 2 J Re r 2 Re 1 dr 


Substitution of (13) into (14) yields 


c 2 * 

rcj = -tt — 
cr * 


1 d^_ dhp^ 

~rlr + dz 2 ' 


(13) 

(14) 

(15) 

(16) 


equations governing the slreamfuncon .^S ? jjSir®traSa1i5Sha». since rhe 
procedure 12 as follows. 

1. Guess the free surface shape for the initial ^ ‘terate ^ transfor ming the governing 

2. Obtain the approximate temperature and y lidna) y domain via a non-orthogonal 

equations and boundary conditions to a method, 

transformation and solve them using *P S transformed momentum equation. 

;s:— 

5 Return totp" 2. Repeat until convergence is obtained by satisfing all equations and 
boundary conditions to a specified degree of accuracy. 

The details of this numerical procedure are discussed below. 

3 . 1 . Non-orthogonal transformation 

(17) 


n = :. 


{-I 


■R(z. t) 

For each iteration time r this transforms the domain 

— A/2 <:< A/2, 0 < r < R(z. t) 
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onto the rectangle 


It then follows that 


— A/2 <2 < a/2, 0<£<I. 


d_ 

dr 


l_d_ 

R dC 


dz 


( 18 ) 


o T ,di»«'Sl f « "" ' m '*7*° » ilh • »- 

^ — »£* 8 o.^ixr;7 e iT™ d " nn8 ,te o " ,cr “ m, ' OT * 

_}_fcifldw Clj/ doj\ uw 1 / d 2 (n r 2 ,„ ■) 32 \ 

1 / Z.l. ^ J -V . _ 


£jA i /d 2 7* 


d z T 

+ A^r + B~ + C^ 


l^C^L) 

c C crjd£ J ’ 


where 


A,(i™ y + _l 

U a?y « ji 

L \RfyJ r c > rj 


(20) 

(21) 




B* = s — £_ t 

R 2 C 

r 2 CcR 

r d n ' 

The boundary and symmetry conditions (5)-(I0) are now: 
at 7 = ± A/2, 


' = 0, «,,*£* 


T=0; 


at { = 0, 


at C= 1, 


s ftf*’ 

$=0, 02 = 0, dT/dZ = 0; 

ip = 0, 


+ f,_ i«v» r mvv<‘iT 

«Un) « a, ^ + L 2 ”U,J J* + *UJ « + L 1+ (^J J v 

1 fi £T c '* f^T lt'Bctni 

(1+(<3B/C27) 2 ) 1 ' j [b ^ a, R dtJ -^JJ + B i {T-r x ) = 0, 


( 22 ) 

(23) 

(24) 

(25) 

(26) 
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G „ , : 2Re ~' 

r i cu / cr y 

( dw 

1 cR £w\ ( 

1 + ((<5K/c>»jJ 2 

lRX + \dn) 1 


R drj dZ ) i 


= Re~ l (Co l — T) 


[1 +(cR/d>j) 2 J 3 ' 2 


( i+mdnl d*R\ 

V R cr, 2 )' (27) 


Equations (24) and (25) are the kinematic condition for the surface (i.e. the condition that a point 
on the surface remain on the surface) and the balance of tangential forces respectively. The 
complete streamfunction, velocity and temperature fields can be determined for a given axisym- 
metnc surface shape (e.g. starting with an initial guess) using only ( !9)-(26). Then, as long as the 
balance-of-force condition (27) is not well satisfied, the force balance residual provides a basis for 
determining an improved estimate of the free surface shape. This procedure continues until the 
convergence criteria (to be defined later) are satisfied. 


3.2. Solution method 

iq 2i 

Various methods can be devised to obtain a steady solution of equations For the sake 

of simplicity a pseudo-unsteady method in association with a semi-implicit time discret iza tion is 
used. 

We first consider the following system of equations: 


‘I + ( t '-SS vl ) I ' +s - 0 ’ 

CCD ( 1 \ 

r +s, *°- 

o. 

where 


(28) 

(29) 

(30) 


. d cij/ 8 \ 

L 1 d c \ u 

R 2 Z \drf c( dtj) R£' 




c n 


S <? 

& ^ 

V 2 =~ — T? 

c n oC 2 \cC/ 


\la ctjcC 


S =-(—-C — ') — 1 — 

2 Re\R 2 C 8tjdc) + Re 2 R dC 


S 3 = 



+ «Cw. 


(31) 

(32) 

(33) 

(34) 

(35) 

(36) 


We then proceed to solve this system as the (pseudo) time derivatives of T, \j) and to-»0. For 
clarity we present the solution method only for the vorticity (o») transport equation. The 
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discretization in time employs an explicit Adams- Bashforth scheme for the non-linear convective 
terms and the implicit scheme for the viscous terms. Other terms in the equation are treated 
explicitly. For spatial discretization, central differences are used. The ADI form for the vorticity 
equation is then 


2 

At 





l 

Re 



C Z OJ _CO )\ H * 112 

A T^y *+■ R — } 
Cijij 


-hs 2 )">— 0, 


(37) 


^ « 1 - o >?/ 1 2 ) + i[L z a < - i(L 2 ajfff 


l / c 2 aj \ n * l 
Re\crj z Jij 


1 

Re 


<? 2 co ccuV + 1/2 

C\/.j 


+ (S 2 )7)=0. 


(38) 


Here the superscript n + 1/2 denotes the intermediate step associated with the ADI method. 16 The 
velocities of (25) (at the mh step) are taken from the values at the previous step. Thus, except for 
the thermal condition at the free boundary, these are Dirichlet boundary conditions. The heat 
transfer equation is also discretized in an ADI form. The resulting system of discretized equations 
is solved using a factorization method. 16 

We define a steady state to occur when the residuals (dcj/dx, dTjcx , cifr/dx) of the vorticity, 
energy and streamfunction equations are less than 10" 7 , i.e. 


I 

Here F represents the vorticity, temperature or streamfunction, the subscript refers to the spatial 
location and the superscript refers to the iterative step. The numerical solution is second- 
order-accurate in space. 

We have verified this solution method for steady 2D axisymmetric incompressible flow. Table I 
shows four examples of a comparison between the results of our code and a finite element code 
(FIDAP); two for buoyancy- and two for thermocapillary-driven flow. Figure 2 shows a com- 
parison of the surface velocities using our method and the finite element method. One sees that 
the numerical method described here provides a reasonably accurate algorithm for computing 
steady 2D axisymmetric incompressible flows in rigid cavities and for mixed buoyancy- and 
surface-driven flows with free boundaries for the range of surface Reynolds numbers examined, 
0<Re<21 740. 

Having computed the vorticity and streamfunction for a given surface shape, it remains to 
iterate on the condition for the force balance normal to the surface in order to obtain the final 
steady surface shape. In addition, the shape must satisfy the volume constraint (11) and boundary 
conditions (12). 


If w + i _f ", 
i r >j r u 


At 


<1(T 


(39) 


Table I. Comparison of finite difference method with results obtained using finite element code FIDAP 


A 

Pr 

Ra 

Re 

This method 

FIDAP 
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2 
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0 
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6-64 x 10* 2 
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6-64 x 10* 2 

2 

0-0127 

0 

3150 

3 39 

1-7 

3-49 

1*82 
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0 

5905 
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Figure 2. Comparison of surface velocities calculated^ Re - 3 150, Gr = 0 and Pr - 001 27 ustng FIDAT 

Two iterative schemes lor determining the and the pressure 
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Scheme 1 is based on the following principle, ^sha^js ^signed f? 

calculated pressure, velocity an ^ directly from (27)using finite differences in con- 

made. The new surface shape js ^ ! vduau* to check whether the 

junction with Newton's .method. The " e graUl 1) ^ itcration is made using a 

procedure u> 'calculate the following improved estimate of * 
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The above procedure for fields are 
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the direction of the force. The magnitude of the iocai displacement is proportional to this force. 
The surface shape at each iteration is thus modified so as to reduce the residual until condition 
(27) is met. It follows that at each iteration the new position of the surface is given by 

R7* l = /?7 + a£x ; , (44) 

where Ex } is the residual of the force balance equation at the jth surface location and the constant 
coefficient a is determined by numerical experiment. In order to ensure convergence, a should be 
small. If a is chosen to be too small, the convergence is slow and the amount of CPU time used 
increases substantially. If x is too large, the solution will diverge. We found that the values of 
x which led to rapid convergence depend on the product of Re~ l and Cq 1 (see Table II)* 

The change in volume between the mth and (m+ l)th iteration can be found from the volume 
constraint (11) and equation (44) and neglecting higher-order terms, i.e. 

*.M2 

ExjR^dn^O. (45) 

tf _ A/2 

The pressure constant /. is contained in Exj and is obtained by satisfying (45). Even then the liquid 
bridge may still change the volume slightly at each iteration owing to numerical error and 
higher-order effects. These small changes can accumulate and eventually result in a gross error. 
To prevent this, formula (44) is modified to 

(v y/* 

*7 +l =*7(“) +*£*/• (46) 

3.3 . Convergence behaviour 

The rate of convergence of the two methods is shown in Figure 3 for Re — 695 and 2082. In both 
cases the shapes of the curves clearly reflect the procedure used. The curves show two distinct 
segments. The abrupt rise in the residual in the late stages of the calculation is caused by requiring 
a more accurate solution (smaller residual) for the vorticity and streamfunction. 

For the calculations carried out here, the first outer iteration scheme converged slightly faster 
than the second. This is shown in Table III. If the guess for the initial iterate is good, the 
computation times are about equal. For a poor guess, however, the first scheme is 1-5 times faster. 
While we were able to obtain convergence for all cases attempted using scheme II, this was not the 
case for scheme I. We found the first scheme to be quite sensitive to the initial choice of the 
pressure constant a. If the value of k is 'physically unreasonable’, the solution diverges. In order to 
take advantage of the speed afforded by scheme I, we employed the second scheme to calculate 
the value of /. for the first few iterations. This value is then taken as an initial guess for the first 
scheme, which is used for the remaining outer iterations. 

The convergence of the solution was also checked by varying the spatial resolution of the mesh. 
This was particularly important at higher values of Re where, owing to the space-centred 
differences, the streamfunction was prone to exhibit ‘wiggles’ if the grid Peclet number 
Pe %Tld = U *A/k (where A is the interval between two grid points) exceeded 2 in the vicinity of the 


Table II. Optimal values of ct as a function of Re l C 0 1 
Re~ l Co' HO 2 0-2 0*1 006 

i 10' 5 !0“ 4 10~ J 2 x 10~ 3 3-3 x I0" 3 
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ITERATION NUMBER 

(a) 

Figure 3. Comparison of convergence 



ITERATION NUMBER 


(b) 

for (a) scheme I and (b) scheme II 


Table III. Sample comparison of CPU times for schemes I and II and 
for the type of initial guess. The parameters are A = 2, Pr = O023, 
Mu = 48, Gr = 76, C 0 = 4-5 x 10 _i and Bi- 100. A ‘good guess’ refers 
to the equilibrium shape of the bridge and a ‘poor guess’ specifies the 
surface shape as R(;)= i + 0006 sin (>(r + 1)]. 



Computation time (s) 

Initial guess 

Scheme I 

Scheme II 

Good 

70 

72 

Poor 

130 

332 


disks. The wiggles are present for V r x JV- = 26 x 61 and are not eliminated until Af,= 10L As 
observed by Ryskin and Leal, 15 an increase in the mesh resolution was found toeliminate this 
problem. We attempted to eliminate the wiggles and avoid the need for mesh sfinement by 
employing second- and third-order upwind schemes for convective terms. We fowl that while 
the wiggles were certainly eliminated for iV r x N s = 26 x 61 (see Figure 4), the meshttiil needed to 
be refined in order to obtain grid convergence. Since the end result was the same, we concluded 
that the centred difference scheme was preferable. A mesh of N r x N s = 26 x 51 was found to be 
sufficient for the results presented here with Reynolds numbers in the range 0 < iSe< 10000. 

We have used the method described earlier to examine the influence of various parameters on 
momentum and heat transport and meniscus shape. In addition, a favourable comparison for full 
and half-zones has been made between the results obtained using the combined scheme and those 
of Hyer et ai 11 who employed a finite element scheme. 
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F.gure 4. Comparison of results obtained using (a) N.-6I, central finite differences, |b) N, -61. third-order upwinding. 
(cl iV, * lot, central finite differences and Id) ,V.= 101, third-order upwinding 


The effects of varying the temperature difference A T (keeping Gr/Re fixed), the ReynoidT 

number Re, the Biot number Bi, the aspect ratio A and the Grashof number Gr are presented in 
the following section. 


4. RESULTS 

4 1. Effects of an increase in AT at fixed Gr/Re 

Figure 5 shows the meniscus shape, dimensionless streamlines and isotherms calculated for 
a fluid with the properties of Gd 5 GajO l2 , Pr = 4-67, Bi=l00, Re = 2 1*4, 107 215 and 1073 
(A/a = 100, 500, 1002 and 5010) and Gr/Re-Q- 025. In all cases the surface tension decreases with 
increasing temperature. This creates a lorce tangent to the surface which drives the melt towards 
the disks and results in the formation of two toroidal rolls with opposite senses of rotation. The 
gravitational acceleration results in a surface shape that bulges out below z = 0 and necks in 

above it. This causes an asymmetry in the structure of the two rolls, with the lower roll being more 
intense. 

At R<? = 21*4 (Ma = 100) heat transport is governed by conduction. In the centre, heat is 
transported through and along the surface {Bi = 100), which results in a significant component of 
the temperature gradient perpendicular to the surface. In the vicinity of the disks the transport of 
heat is mostly through conduction along the surface. At higher values of Re convection is 
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intensified. The flow velocities reach the,. »!»« i * “S^l^rTflow 

-entral Doruon of the surface and the disks. In the central region of the melt the outward no 

towards the surface bnngs cooler melt from the interior and causes a steep temperature grad.en 

perpendicular to the surface. Towards the end walls the frequency 

the surface increases This indicates a steeper temperature gradient parallel to the surtace. 

ngtc 6 shows ,h. effect of increasing Cr and Re a. fined Gr/Re |.M» "‘ h ^ 

properties of molten silicon. As AT is increased, the increase m bao»anc,^n«n do. m he upper 
region results in a larger upper cell which extends into the lower half. Note that the ratio ol he 
magnitudes of K,. to L* '"creases as Gr and Re are increased. It is £ 

GriRe is small, if the ratio remains constant but the magnitudes of Grand Re are ^rease^t te 
buoyanev-driven effect manifests itself more at higher Gr- an e-va u . which acts 

buoyancy effect due to the (vertically) thermally unstable s.tuanon^ 
together with the radial temperature gradients throughout the hqutd bndge u !£*•«» 
ward motion of colder fluid in the region midway between the axis and l the i surf . u^ 

buoyancy-driven flow becomes pervasive throughout the system and confines the more intense 
(thermocapillary-driven) cell to a small region in the lower half. 

The effect of increasing AT on the free surface shape is insignificant for these cases. 
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Figure 6. Meniscus shape, streamlines and isotherms calculated for a fluid with the properties of molten silicon, 
Pr = 0-023, Bt = 100 and Gr 0-036: U)*<f=I39. c„ = 23xlO~ J ; <b)tf*=l390, c, = I 4 x 10~ 3 ; {c)/*c = 2082, 

fy- I I x 10" 3 ; c t -0 09 for all three cases 


4.2. Effect of Biot number 

The effect of Bi on the liquid bridge shape is slight for the cases examined. The isotherm 
distribution is modified, however, when Bi is increased from 10 to 100 (see Figures 5(b) and 7), In 
the central part of the surface the temperature gradient is almost perpendicular to the surface in 
both cases. For Bi = 100 the gradient is steeper and the temperature in the central region is higher, 
which results in a steeper temperature gradient parallel to the surface near the disks. There is 
a slight decrease in flow intensify for the Bi= 10 case. 
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Figure 7. Effect of Bi (compare with Figure 5) for Si - 10, Pr- 4-667, Re l ^7, C 0086 resw^eiv ^ ^ 
8 streamfunciion and temperature contour intervals are 19 x 10 and 0-086 respectively 


= 2-7. The 


4.3. Effect of aspect ratio (A) 

For otherwise identical physical conditions a change in aspect ratio of the liquid bridge has 
a significant effect on the meniscus shape. This is illustrated in Figure 8 by a comparison of liquid 
bridges with A = 20 and 3. The shapes are qualitatively similar but the amph ude of he 
deformation from a right circular cylinder is much larger for the higher aspect ratio. For the 
large-aspect-ratio bridge there is a relative decrease in flow intensity in the upper P art as thc flow 
conforms to the increased curvature of the surface of the longer liquid bridge. The flow intensity 
of the lower cell is greater for the longer aspect ratios and the isotherms exhibit correspondingly 
more distortion. Results obtained for A = £ confirm this trend (note that with a Bond number 
Bo = pgR 1 h o = 0 - ? the A = 3^2 bridge is close to the stability limit ). 


4.4 . Pure thermocapillary flow 

In the absence of gravity Gr = B o = 0 and the flow is driven only by surface tension gradients. 
Figure 9 shows the isotherms and streamlines computed for a fluid wlt ^ the P^ 0 ^* eS of 
silfcon for values of Re up to 21 740. We computed cases with Re= 1390, 4350, 8695 and 21 740. 
At Re 5 4350 we found secondary cells. There is an increase in secondary flow intensiy at higher 
values of Re (see Figures 9(b) and 9(c)). Even at the higher values of Re the isotherms are only 

^'FoMhe'finilfgravity case (Gr = 76) the effects of surface shape and buoyancy-driven flow are 
manifested (see Figure 10). Only two cells are evident in the finite gravity case -The larger ceU 
appears to be a combinat.on of a primary thermocapillary cell and a downward flow due to 
buoyancy near the axis, which extends into the central region of the lower half and interacts with 
a secondary flow cells associated with the thermocapillary flow. The primary thermocapillary- 
driven cell in the lower half is considerably smaller in extent than it is for the zero-gravity case. 
A comparison between zero-gravity and finite gravity conditions illustrates the interaction 
between the weaker buoyancy-driven flow, the secondary thermocapillary cells (see Figures 9 
and 10) and the surface shape. 
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T ^-0 

(a) 


WO 



Figure 8. Effect of aspect ratio: (at A = 2-5; (b) A = 3; Bi= 100, B o = 01, Fr = 4 667, Re= 107, C~ 1 = 70 C = 0468 Gr = 2-7 
rhe sireamfunction contour intervals are 1 2x KT 3 and 2-3 x KT 3 for (a) and (b) respectively and the temperature 

contour interval is 0-087 for both cases 


5. CONCLUSIONS 

Steady solutions to the free boundary problem for a non-isothermal liquid badge have been 
obtained using a Picard-type iterative scheme with a semi-implicit space-centred finite difference 
scheme. The method was applied to a variety of problems with Reynolds numbers in the range 
0< Re <22000 and for Pradtl numbers of 0 023 and 4-67. The method was found to perform well 
for this range of parameters examined and compared well with results obtained for buoyancy- 
dnven and surface-tension-driven Row using a finite element method. At high Re we found that 
the solution was sensitive to the spatial resolution and that a larger number of grid points were 
necessary in order to avoid ‘wiggles’ in the solution. For Re< 10000 we found N r x jV. = 26 x 51 
points to be sufficient, but required up to 101 points in the ^-direction for Re = 21 740. While it 
was possible to remove the wiggles using upwind differences, the mesh refinement was still 
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Pr-^7? i,n 7^ llnes and ,solhern ’ s “Iculaled for a fluid with the properties of mollen silicon, 

" 0023, a -100, Gr-76 and Re = 8695. The streamfunction contour interval is 46x I(T* and the temperature 

contour interval is 0*087 


ATis increased. As expected from the results of previous work, the effect of changing the aspect 
ratio was also seen to significantly affect the flow behaviour. For non-zero gravity the increase in 
convexity of the lower portion of the larger-aspect-ratio bridge has a pronounced effect on the 
flow pattern, which is more intense for the longer liquid bridge. This increase in flow intensity 
appears to occur as a result of a decrease in the form drag of the surface. A comparison between 
the difference in flow intensity between the upper and lower halves reveals that the difference 
increases with increasing A. The form drag in the upper half increases as the magnitude of the 
negative curvature increases. This sfrog s the flow in the upper half relative to the lower half of the 
bridge. 
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Abstract _ 

Liquid bridges appear in a variety of industrial processes, forexamplein the well-known floating-zone 
crystal growth technique. This crystal growth method has received much attention in recent years. 
In particular, there have been a variety of experiments on spacelab missions. These experiments are 
motivated by the fact that the microgravity environment affords the possibility of an increase in the 
stability of the melt meniscus and a reduction in buoyancy-driven convection. However, within the 
spacecraft there is a residual acceleration with variable magnitude and orientation. Under certain 
conditions, the response of the free surface of a liquid bridge to time-dependent residual accelerations 
will lead to zone breakage. In this paper the steady and unsteady behavior of isothermal and non- 
isothermal liquid bridge systems under normal and low gravity conditions is examined. The full non- 
linear governing equations are recast in terms of a stream-function vorticity formulation together with 
a non-orthogonal coordinate transformation. The latter allows an irregular free boundary to coincide 
with a coordinate line (or surface) without the need to solve a coupled set of Laplace equations. The 
resulting equations are discretized using a centered finite difference scheme for space, and an Adams- 
Bashforth-Crank-Nicolson scheme is used for time. The equations are solved by the A.D.I. method 
and a Picard type iteration is used on the boundary condition for the balance of force normal to the 
free surface. For non-isothermal bridges, residual acceleration affects the system by causing internal 
buoyancy flows and fluctuations in the shape of the bridge which interact with the thermocapillary 
flow caused by surface tension gradients. For the cases examined, the shape of the bndge is found 
to be more sensitive to typical spacecraft accelerations than the buoyancy driven flow. The effect of 
thermocapillary flow on the surface shape is found to be small for the range of capillary and Reynolds 
numbers considered. 

1. Introduction 

Liquid bridges appear in a variety of industrial processes, for example the well-known 
floating-zone crystal growth technique [1]. This crystal growth method has received much attention 
in recent years [2- 1 0] . In particular, there have been several related experiments on spacelab missions 
[11-14]. These experiments are motivated by the fact that the microgravity environment affords the 
possibility of an increase in the stability of the melt meniscus and a reduction in buoyancy-driven 
convection. However, within the spacecraft there is a residual acceleration with variable magnitude 
and orientation. Under certain conditions, the response of the free surface of an isothermal liquid bndge 
to time-dependent residual accelerations will lead to zone breakage [13-15]. In this paper we examine 
the interaction between convection caused by the response of the free surface of the zone to oscillatory 
axial residual acceleration and convection due to thermocapillary and internal buoyancy forces. 
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Fig. 1 . The model liquid bridge. 


t mn =o 



Fig. 2. The initial steady dimensionless stream- 
function and temperature corresponding to Re = 2899, 
Gr = 0.98, A = 4 and Pr = 0.0127. 


2. Formulation 

2.1 Description of the Model 

A cylindrical liquid zone (see Fig. 1) is contained between two parallel coaxial circular rigid 
disks (radius =R q) separated by a distance L The liquid is a non-isothermal Newtonian fluid for which 
the Boussinesq approximation holds. The bridge is held between the disks by surface tension. The 
free surface of the bridge is a gas-liquid interface and is described by r = R(z,t). Each disk is maintained 
at a constant temperature To. Surface heating is provided through a parabolic function, T„, which is 
a function of the axial coordinate and is an approximation to the heating profile associated with typical 
floating zone crystal growth experiments. The heat transfer coefficient at the free surface is denoted 
by h. In addition, we make the assumptions that the residual acceleration is parallel to the cylinder 
axis, the velocity and temperature field and the deformation of free surface are axisymmetric, and we 
take the surface tension at the free surface to be a linear function of the temperature. Motion of the 
end disks perpendicular to the axis may also occur in practice. The restriction to axial acceleration 
precludes an analysis of the effects of such motions. 

The governing equations are made dimensionless by scaling length, time and velocity with 
Ro, Ro/U* and U*, respectively. Here U* is a characteristic velocity given by 



M- 


where AT = Tmax - Tmin represents the maximum temperature difference along the surface, [yj is 
the absolute value of the derivative of the surface tension with respect to temperature, and |i is the 



dynamic viscosity. The difference AT is used to non-dimensionalize temperature. 

The temperature maximum in a floating zone occurs between the two ends of the zone. We 
shall take the dimensionless ambient temperature Tmax to be T„,(0), and T min to beT«,(±A/2), where 
A = L/R 0 is the aspect ratio. 

The non-dimensional pressure is 

P* + pog*(t) z 


P = - 


poU 


■*2 


Ro, 


where p* is the dimensional pressure, z is the dimensionless axial coordinate, po is the density cor- 
responding to the reference temperature and g*(t) = go + gisin(2jt/t) is the residual gravitational 
acceleration. 


2.2 Basic Equations 

With the scales presented in 2.1 the governing dimensionless equations in a cylindrical 
coordinate system can be written as 
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where g(t) = g*(t)/go is a time-dependent dimensionless residual acceleration, and 


Re = 


_R 0 U 


„ |T t I ATR o ^ IgolPATR^ 

Ma = — u , Or = 

|IK 


V |IK V 

are, respectively, the Reynolds number, Marangoni number and Grashof number. Here, v is the ki- 
nematic viscosity, k is the thermal diffusivity, (3 is the volume thermal expansion and coefficient. 
The boundary conditions at the rigid end disks are 


u = T = 0, w = w ± (t), atz = ±A-, 


(5) 


where w— (t) is zero if the two disks vibrate in phase (this will admit only odd mode deformations of 
the zone surface [15] ). The symmetry conditions at the centerline r = 0 are 

3w = 3T ( 6 ) 

3r 3r 

The boundary conditions at the free surface r = R(z) take the form 




Fig. 3. The instantaneous dimensionless stream-functions with Re = 2899, Gr = 0.98, A = 4 
and Pr = 0.0127, at (a) 0.52 s, (b) 0.59 s, (c) 1 s, (d) 1.5 s (e) 156 s, (0 2 s, after application of 
an additional 2.5xl0" 2 g, 0.5 Hz, axial acceleration. 


g(t)z , 2Re 4 r d u t pR| 2 8w dR /8w ( 9u \ + 
F 2 . /3R j 2 .dr \3zj dz dz l Or dz ) 

\5z / 


Re' 1 (Cp 1 - t) 




’i.f^| 2 l[^ + §w| _[ /or_\ 2 | z far aR§T\ 

\dz / \5z dr I dz \3r dz ) |dz ) \3z dz dr j ’ 



( 7 ) 

( 8 ) 

( 9 ) 

( 10 ) 




and Yo is the mean surface tension and k is the thermal conductivity. The force balance conditions 
normal and tangent to the free surface are given by eqs. (7) and (8) respectively. Equation (9) is the 
kinematic boundary condition at the liquid-gas interface. The thermal boundary condition at the 
interfaceis given by equation ( 1 0) in which the equivalent heat transfer coefficient, h, models the effect 
of the radiant and convective heat transfer between the bridge and the surrounding environment. The 
constant X in (7) represents a dimensionless reference pressure difference across the interface which 
for this system is determined by the following constant volume constraint [9,16] 

fj 

V = I 7 tR 2 (z)dz = V 0 = constant. (H) 

2 


Finally, the condition that the contact lines between the liquid end disks are fixed is 

R = 1 at z = ± A/2. 



Fig. 4. The dimensionless temperature field with Re = 2899, Gr = 0.98, A - 4 and 
Pr = 0.0127, at (a) 0.52 s, (b) 0.59 s, (c) 1 s, (d) 1.5 s (e) 1.56 s, (0 2 s, after 
application of an additional 2.5x1 0"^ g, 0.5 Hz, axial acceleration. 




2.3 Solution Method 

In the present investigation, the governing equations are recast in terms of a stream-function 
vorticity formulation. The stream-function is defined by 


u = 


1 9Xp 
r dz 
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l dV 
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A non-orthogonal coordinate transformation, 

Z 




(14) 


R(z,t) ’ 

allows an irregular fiee boundary to coincide with a cylindrical coordinate line (or surface) without the 
need to solve a coupled set of Laplace equations [17,18], The resulting equations are discretized 
following a semi-implicit difference scheme and solved by the A.D.I. method. The conditions for force 
balance tangent to the surface and kinematic condition at the free surface are solved along with the 
Navier-Stokes and continuity equations. The condition for the force balance normal to the surface is 
used together with an “outer” iterative procedure to determine the free surface shape. 

The unsteady free boundary problem for a cylindrical liquid zone is solved as follows. The 
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Fig. 5. The instantaneous dimensionless stream- 
function, with Re = 0, Gr = 0, A = 4 and Pr = 0.0127, 
at (a) 0.59 s, (b) 1 s, (d) 1.5 s (e) 2 s, after application 
of an additional 2.5x10" 2 g, 0.5 _Hz, axial acceleration. 
Note that for this case only, v P= T/ kRo , where'P is the 
dimensional stream-function cm 3 s' 1 . 



Fig. 6. Tolerable acceleration [g] vs. Frequency [Hz] 
predicted by the ID model [15] for (1) breakage and 
(2) 10% shape change, for an isothermal liquid zone 
with A = 4 and the thermophysical properties of 
indium. 



initial conditions correspond to either zero or finite steady residual acceleration situations with a steady 
thermocapillary flow. These states are calculated using a method close to that described below [16]. 
For the unsteady calculation our solution scheme is similar to that used by Kang and Leal [17] and 
Ryskin and Leal [18]. The following Picard iterative procedure [19] is adopted. 

1. guess the free surface shape for the new timestep; 

2. obtain the approximate temperature and stream-function, vorticity and velocity fields by transform- 
ing the governing equations and boundary conditions to a circular cylindrical domain via a non- 
orthogonal transformation and solve them using a semi-implicit method, 

3. obtain the pressure at the free surface by integrating the transformed momentum equation; 

4. use the condition for the balance of force normal to the free surface to decide how to update the free 
surface location; 

5. return to step 2. Repeat until convergence is obtained by satisfying all equations and boundary 
conditions to a specified degree of accuracy for this timestep. 

3. Results and Discussion 

The following results were obtained for a liquid zone corresponding to the physical properties 
of molten indium subject to an axial acceleration with a frequency of 0.5 Hz. Fig. 2 depicts the initial 
state of the system . A steady axial acceleration of magnitude 10 4 g (lfr 3 m s’ 2 ) acts along the negative 
z-direction. Two equidimensional toroidal rolls indicate that surface-driven flow is dominant (Re = 
2899, Gr = 0.98). The isotherm distribution shows that heat transfer is mainly by conduction, although 
some distortion of the isotherms by the flow is evident. For indium, Lind [20] has reported that the 
surface tension increases with increasing temperature i.e. Yr > 0. (Note that surface contamination may 
have affected the measured temperature dependence of surface tension in this case.) Thus, the flow 
direction at the surface is toward the center (i.e. the higher temperature region). Figs. 3 and 4 illustrate 
the effect of an additional acceleration component which varies sinusoidally with a frequency of 0.5 
Hz. Figure 5 depicts the response of the zone to the same disturbance, but with Gr = Ma = 0. (Note 
that, for this case, k/L rather than fyj AT/p was used as the velocity scale.) In another case with Re = 
2899, Gr = 0.98 and the surface constrained to be a circular cylinder, no observable response occurred. 
Clearly the system is more sensitive to the effects of free surface motion than internal buoyancy. Given 
that the value of steady acceleration used is extreme for spacecraft acceleration environments [2 1 ], we 
may conclude that for systems where internal buoyancy-driven effects are not manifested (in this case 
because they are swamped by the surface-driven flow) it suffices to examine the response of the free 
surface only. Furthermore, a comparison of our full axisymmetric results with those obtained with a 
simplified ID isothermal model indicates that, at least for the conditions examined, the ID model may 
be used to reliably predict liquid zone (isothermal and nonisothermal) sensitivity. The ID model is 
described in detail in [ 15]. Fig. 6 shows the sensitivity of an indium liquid zone to axial acceleration. 
The curves are based on results obtained using two sensitivity criteria. The first is determined by 



breakage of the bridge, the second is whenever the bridge shape changes by more than 1 0% of its static 

shape, i.e. R(z,t) - R(z,0) = .lR(z,0). 
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Liquid bridges appear in a variety of industrial processes, for example in the well-known floating-zone 
crystal growth technique. This crystal growth method has received much attention in recent years. 
In particular, there have been a variety of experiments on spacelab missions. These experiments are 
motivated by the fact that the microgravity environment affords the possibility of an increase in the 
stability of the melt meniscus and a reduction in buoyancy-driven convection. However, within the 
spacecraft there is a residual acceleration with variable magnitude and orientation. Under certain 
conditions, the response of the free surface of a liquid bridge to time-dependent residual accelerations 
will lead to zone breakage. In this paper the steady and unsteady behavior of isothermal and non- 
isothermal liquid bridge systems under normal and low gravity conditions is examined. The full non- 
linear governing equations are recast in terms of a stream- function vorticity formulation together with 
a non-orthogonal coordinate transformation. The latter allows an irregular free boundary to coincide 
with a coordinate line (or surface) without the need to solve a coupled set of Laplace equations. The 
resulting equations are discretized using a centered finite difference scheme for space, and an Adams- 
Bashforth-Crank-Nicolson scheme is used for time. The equations are solved by the A.D.I. method 
and a Picard type iteration is used on the boundary condition for the balance of force normal to the 
free surface. For non-isothermal bridges, residual acceleration affects the system by causing internal 
buoyancy flows and fluctuations in the shape of the bridge which interact with the thermocapillary 
flow caused by surface tension gradients. For the cases examined, the shape of the bridge is found 
to be more sensitive to typical spacecraft accelerations than the buoyancy driven flow. The effect of 
thermocapillary flow on the surface shape is found to be small for the range of capillary and Reynolds 

numbers considered. 

1. Introduction 

Liquid bridges appear in a variety of industrial processes, for example the well-known 
floating-zone crystal growth technique [ 1 1 . This crystal growth method has received much attention 
in recent years [2- 1 0] . In particular, there have been several related experiments on spacelab missions 
[11-14], These experiments are motivated by the fact that the microgravity environment affords the 
possibility of an increase in the stability of the melt meniscus and a reduction in buoyancy-driven 
convection. However, within the spacecraft there is a residual acceleration with variable magnitude 
and orientation. Under certain conditions, theresponse of the free surface of an isothermal liquid badge 
to time-dependent residual accelerations will lead to zone breakage [13-15]. In this paper we examine 
the interaction between convection caused by the response of the free surface of the zone to oscillatory 
axial residual acceleration and convection due to thermocapillary and internal buoyancy forces. 
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Fig. 1. The modelliquid bridge. 
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Fig. 2. The initial steady dimensionless stream- 
function and temperature corresponding to Re = 2899, 
Gr = 0.98, A = 4 and Pr = 0.0127. 


2. Formulation 


2.1 Description of the Model 

A cylindrical liquid zone (see Fig. 1) is contained between two parallel coaxial circular rigid 
disks (radius = R 0 ) separated by adistance L The liquidis a non-isothermal Newtonian fluid for which 
the Boussinesq approximation holds. The bridge is held between the disks by surface tension. The 
free surface of the bridge is a gas-liquid interface and is described by r =R(z,t). Each disk is maintained 
at a constant temperature T 0 . Surface heating is provided through a parabolic function, T^ which is 
a function of the axial coordinate and is an approximation to the heating profile associated with typical 
floating zone crystal growth experiments. The heat transfer coefficient at the free surface is denoted 
by h. In addition, we make the assumptions that the residual acceleration is parallel to the cylinder 
axis, the velocity and temperature field and the deformation of free surface are axisymmetric, and we 
take the surface tension at the free surface to be a linear function of the temperature. Motion of the 
end disks perpendicular to the axis may also occur in practice. The restriction to axial acceleration 
precludes an analysis of the effects of such motions. 

The governing equations are made dimensionless by scaling length, time and velocity with 
Ro, Rq/U* and U*, respectively. Here U* is a characteristic velocity given by 


* 

U 



where AT = Tmax - Tmin represents the maximum temperature difference along the surface, fyj is 
the absolute value of the derivative of the surface tension with respect to temperature, and ji is the 



dynamic viscosity. The difference AT is used to non-dimensionalize temperature. 

The temperature maximum in a floating zone occurs between the two ends of the zone. We 
shall take the dimensionless ambient temperature T max to be T«(0), and Tmin to beToo(±A/2), where 
A = L/Ro is the aspect ratio. 

The non-dimensional pressure is 

P* + pog*(0 z D 

P= = K 0’ 

PoU 2 

where p* is the dimensional pressure, z is the dimensionless axial coordinate, po is the density cor- 
responding to the reference temperature and g*(t) = go + gisin(2rc/t) is the residual gravitational 
acceleration. 


2.2 Basic Equations 

With the scales presented in 2.1 the governing dimensionless equations in a cylindrical 
coordinate system can be written as 
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where g(t) = g*(t)/go is a time-dependent dimensionless residual acceleration, and 

Rc= MJl , MaJj&, Gr JSS^L. 

V J1K 

are, respectively, the Reynolds number, Marangoni number and Grashof number. Here, v is the ki 
nematic viscosity, K is the thermal diffusivity, (3 is the volume thermal expansion and coefficient. 
The boundary conditions at the rigid end disks are 


u= T = 0, w = w ± (t), atz=±^, (5) 

where w^t) is zero if the two disks vibrate in phase (this will admit only odd mode deformations of 
the zone surface [ 15] ). The symmetry conditions at the centerline r = 0 are 




The boundary conditions at the free surface r = R(z) take the form 
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Fig. 3. The instantaneous dimensionless stream-functions with Re = 2899, Gr = 0.98, A = 4 
and Pr = 0.0127, at (a) 0.52 s, (b) 0.59 s, (c) 1 s, (d) 1.5 s (e) 1.56 s, (0 2 s, after application of 
an additional 2.5x1 0 -2 g, 0.5 Hz, axial acceleration. 
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where the capillary number, Biot number and Froude number are, respectively, 
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and Yo is the mean surface tension and k is the thermal conductivity. The force balance conditions 
normal and tangent to the free surface are given by eqs. (7) and (8) respectively. Equation (9) is the 
kinematic boundary condition at the liquid-gas interface. The thermal boundary condition at the 
interface is given by equation ( 1 0) i n which the equivalent heat transfer coefficient, h, models the e ffect 
of the radiant and convective heat transfer between the bridge and the surrounding environment. The 
constant X in (7) represents a dimensionless reference pressure difference across the interface which 
for this system is determined by the following constant volume constraint [9,16] 

fj 

V = I TtR 2 (z)dz = V 0 = constant. (11) 

2 

Finally, the condition that the contact lines between the liquid end disks are fixed is 

R = 1 at z = ± A/2. (12) 



Fig. 4. The dimensionless temperature Field with Re = 2899, Gr = 0.98, A - 4 and 
Pr = 0.0127, at (a) 0.52 s, (b) 0.59 s, (c) 1 s, (d) 1.5 s (e) 1.56 s, (f) 2 s, after 
application of an additional 2.5x1 0"^ g, 0.5 Hz, axial acceleration. 




2.3 Solution Method 

In the present investigation, the governing equations are recast in terms of a stream-function 
vorticity formulation. The stream-function is defined by 


= _i^ 

,V r 3r 


A non-orthogonal coordinate transformation, 

n=z ' 5 = r&' (14) 

allows an irregular free boundary to coincide with a cylindrical coordinate line (or surface) without the 
need to solve a coupled set of Laplace equations [17,18]. The resulting equations are discretized 
following a semi-implicit difference scheme and solved by the A.D.I. method. The conditions for force 
balance tangent to the surface and kinematic condition at the free surface are solved along with the 
Navier-Stokes and continuity equations. The condition for the force balance normal to the surface is 
used together with an “outer” iterative procedure to determine the free surface shape. 

The unsteady free boundary problem for a cylindrical liquid zone is solved as follows. The 
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Fig. 5. The instantaneous dimensionless stream- 
function, with Re = 0, Gr = 0, A = 4 and Pr = 0.0127, 
at (a) 0.59 s, (b) 1 s, (d) 1.5 s (e) 2 s, after application 
of an additional 2.5xl0“ 2 g, 0.5 Hz, axial acceleration. 
Note that for this case only, x F=T/td^o , where 1 ? is the 
dimensional stream-function cm 3 s* 1 . 
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Fig. 6. Tolerable acceleration [g] vs. Frequency [Hz] 
predicted by the ID model [15] for (1) breakage and 
(2) 10% shape change, for an isothermal liquid zone 
with A = 4 and the thermophysical properties of 
indium. 



initial conditions correspond toeither zero or finite steady residual acceleration situations with a steady 
thermocapillary flow. These states are calculated using a method close to that described below [16]. 
For the unsteady calculation our solution scheme is similar to that used by Kang and Leal [17] and 
Ryskin and Leal [18]. The following Picard iterative procedure [19] is adopted: 

1. guess the free surface shape for the new timestep; 

2. obtain the approximate temperature and stream-function, vorticity and velocity fields by transform- 
ing the governing equations and boundary conditions to a circular cylindrical domain via a non- 
orthogonal transformation and solve them using a semi-implicit method; 

3. obtain the pressure at the free surface by integrating the transformed momentum equation; 

4. use the condition for the balance of force normal to the free surface to decide how to update the free 
surface location; 

5. return to step 2. Repeat until convergence is obtained by satisfying all equations and boundary 
conditions to a specified degree of accuracy for this timestep. 

3. Results and Discussion 

The following results were obtained for a liquid zone corresponding to the physical properties 
of molten indium subject to an axial acceleration with a frequency of 0.5 Hz. Fig. 2 depicts the initial 
state of the system . A steady axial acceleration of magnitude 10 4 g (10* 3 m s' 2 ) acts along the negative 
z-direction. Two equidimensional toroidal rolls indicate that surface-driven flow is dominant (Re = 
2899, Gr = 0.98). The isotherm distribution shows that heat transfer is mainly by conduction, although 
some distortion of the isotherms by the flow is evident. For indium, Lind [20] has reported that the 
surface tension increases with increasing temperature i.e. Yr>0. (Note that surface contamination may 
have affected the measured temperature dependence of surface tension in this case.) Thus, the flow 
direction at the surface is toward the center (i.e. the higher temperature region). Figs. 3 and 4 illustrate 
the effect of an additional acceleration component which varies sinusoidally with a frequency of 0.5 
Hz. Figure 5 depicts the response of the zone to the same disturbance, but with Gr = Ma = 0. (Note 
that, for this case, k/L rather than fyj AT/(i was used as the velocity scale.) In another case with Re = 
2899, Gr = 0.98 and the surface constrained to be a circular cylinder, no observable response occurred. 
Clearly the system is more sensitive to the effects of free surface motion than internal buoyancy. Given 
that the value of steady acceleration used is extreme for spacecraft acceleration environments [2 1 ] , we 
may conclude that for systems where internal buoyancy-driven effects are not manifested (in this case 
because they are swamped by the surface-driven flow) it suffices to examine the response of the free 
surface only. Furthermore, a comparison of our full axisymmetric results with those obtained with a 
simplified ID isothermal model indicates that, at least for the conditions examined, the ID model may 
be used to reliably predict liquid zone (isothermal and nonisothermal) sensitivity. The ID model is 
described in detail in [15]. Fig. 6 shows the sensitivity of an indium liquid zone to axial acceleration. 
The curves are based on results obtained using two sensitivity criteria. The first is determined by 


breakage of the bridge, the second is whenever the bridge shape changes by more than 1 0% of its static 

shape, i.e. R(z,t) - R(z,0) = .lR(z,0). 
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